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ABSTRACT 

An analytical approach to the solution of the whipping response of a submerged submarine 
to the pulsation of a nearby explosion bubble is presented Particular aspects of energy 
dissipation (damping) are addressed in regards to their effects on the overall loading forces 
and dynamic response of the submarine. Particular damping mechanisms discussed 
include hull damping (hull material and structural), internal damping (internal structures 
and sloshing liquids), and external damping (hydrodynamic damping including wave 
radiation and viscous fluid effects). 

A flexible analytic model is developed, using the fundamentals of vibration and hydro- 
mechanics. A finite-element beam model is used to represent the flexural structure of the 
hull. Internal structures are modeled as separate dynamic systems, using both discrete and 
modal superposition approaches. Onboard liquids are modeled using a mechanical- 
equivalent system based upon a potential flow solution of the liquid free-surface. External 
hydrodynamic forces are modeled using a modified Morison formulation, with fluid 
velocities and accelerations calculated based upon a popular explosion bubble model. The 
equations of motion for the composite dynamic systems are solved in the time-domain 
using a modified Newmark direct time integration scheme, with iterations at each time 
step accomplished using a modified Newton-Raphson method 

The analytic model is implemented on a personal computer, and used to numerically 
investigate various damping mechanisms associated with dissipation of whipping energy 
Additionally, the analvtic model is compared to explosion-induced whipping tests 
conducted on the U.S. Navy test platform "Red Snapper". The analytic model compared 
well to the experimental data for early-time response (through the 2nd bubble period), but 
tended to over predict response for later times. 

Potential weaknesses with the analytic model are discussed in light of the comparisons 
with experimental data, and recommendations for future analytic work in the area are 
made. 

Thesis Supervisor: Professor Alan Brown, Department of Ocean Engineering 
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1. INTRODUCTION 



1.1 Underwater explosions and ship whipping. 

Designers of naval ships have long been concerned with the effects of underwater 
explosions. Historically, most attention has been paid to the effects of contact charges and 
charges of near to intermediate standoff - . Protection methods to these effects have 
typically included heavy side-protective systems to minimize possibility of hull rupture and 
shock-mounting of vital machinery and equipment to minimize damage caused by the high 
frequency shock wave resulting from the explosion. 

While hull rupture and internal shock wave effects are still considered of major 
importance to ship survivability, increasing attention has been paid in recent years to the 
whipping of ship hulls caused by charges of near to intermediate standoff Whipping can 
be defined as the transient, flexing motion of the whole hull of a ship in its vibrational 
modes corresponding to the lowest natural frequencies 1 . Such motion, if of sufficient 
magnitude, can cause hull girder buckling, tearing, or other loss of hull girder strength 

Numerous examples of catastrophic or near-catastrophic whipping can be cited 
throughout naval history. As recently as 1991, the U S. Navy Guided Missile Cruiser 
USS Princeton (CG-59) struck a floating contact mine in the Persian Gulf Although the 
local hull rupture caused by the mine explosion (near the bow) was not a threat to the 
survival of the ship, the subsequent whipping of the hull caused near-catastrophic hull 
girder damage near the stern of the ship 

When an underwater explosive detonates, the solid explosive material suddenly 
reacts, leaving behind gaseous products at very high temperature and pressure Almost 
immediately, the shock wave, produced due to the sudden discontinuity in pressure, 
propagates radially, approximately with the speed of sound in water. This shock wave 
carries with it approximately one-third of the energy released during the reaction 2 . It is 
this shock wave that is typically responsible for localized hull rupture and damage to 
internal equipment. 

Left behind as the shock wave propagates through the fluid medium is a cavity of 
gaseous reaction products at high pressure This cavity subsequently expands (as a 
bubble) to relieve the difference in pressure, accelerating the surrounding fluid. The 
bubble continues to expand beyond the point of hydrostatic equilibrium (due to the inertia 
of the surrounding fluid) until a point of dynamic equilibrium is reached The bubble then 

'Hicks. A.N.. "Explosion Induced Hull Whipping". . Usances in Marine Structures. Admiralty Research 
Establishment. Dunfermline. Scotland. UK. 1986, p. 390. 

2 Ibul. p. 392. 
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reverses, continuing to contract until dynamic equilibrium is again reached, where it 
quickly rebounds and again begins to expand. This bubble expansion and contraction 
continues until the energy of the reaction is fully dissipated. As the bubble rebounds, it 
greatly accelerates the surrounding water, generating a substantial pressure pulse 
(commonly known as the bubble pulse) This bubble pulse can impart significant loads on 
structures in the vicinity. 

Thus, for certain underwater explosions, a ship can experience loading, not only 
from the radiating shock wave, but also from the quasi-periodic bubble pulse. Although 
the whipping motion of a ship hull is characterized by free vibration of the hull in the 
water, it is this quasi-periodic loading which can reinforce the free vibration and magnify 
the response of the hull to whipping, particularly when the period of the bubble pulse is in 
the vicinity of the natural periods of hull flexure. 

1.2 The importance of damping in ship whipping. 

The first, most easily observed effect of damping on structural vibration is in the 
decrease in amplitude during free vibration. Figure 1-1 illustrates the effect of damping 
on the free vibration of a cantilever beam. If the cantilever with some damping is 
deformed and then released, it will oscillate with a regular period, but the amplitude of the 
oscillation will decrease with time. Without damping, the cantilever would continue to 
vibrate without the decrease in amplitude. In fact, the free vibrations of damped structures 
will normally be reduced at a rate that may be used as a measure of the amount of 
damping in the system A well known measure of damping, known as the logarithmic 
decrement, 5, is related to the ratio of the /rth to the (n^Njth cycle amplitudes by: 

s = -U(-^-) (i-D 

N a „.n 

where A n is the amplitude of the /rth cycle and A n+N is the amplitude of the (//+A0th 
cycle. Similarly to a cantilever, the amplitude of free vibration of a whipping ship will 
decrease with time due to the presence of damping. 

A second, and perhaps more important effect of damping is on the steady state 
amplitude attained when a structure is excited by a harmonically oscillating force. As 
discussed in any text on mechanical vibrations, the maximum magnitude of the response of 
any structure subjected to harmonic excitation (or periodic excitation that can be 
represented as a Fourier series of harmonic terms) is determined by the relation of the 
frequency of excitation to the natural frequencies of the structure, as well as the stiffness. 
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Figure 1-1 Effect of damping (with time) on the free vibration 
of a cantilever. 



damping and mass of the structure (and, of course, the magnitude of the excitation). As 
shown in figure 1-2 (illustrated for a single degree of freedom system subjected to 
harmonic excitation), the response is predominantly controlled by stiffness when the 
excitation frequency (co) is much less than the natural frequency (co n ) of the structure 
(elastic forces dominate the equations of motion due to the small acceleration and inertia 
forces). When the excitation frequency is much greater than the natural frequency, the 
response is predominantly controlled by mass (inertia forces dominate the equations of 
motion) When the excitation frequency is close to the natural frequency (near 
resonance), however, the elastic and inertia forces approximately balance, and the major 
demand on the excitation forces is to overcome system damping, implying that the amount 
of system damping controls the dynamic response. In this region, when damping is low, 
the response magnitude is high, and when damping is high, the response amplitude is low. 
Thus, the response near resonance can be very much greater than the static response (co 
approaching zero), particularly if the damping is very low 

As discussed previously, the bubble pulse can generate low frequency, quasi- 
periodic loading on the ship. For certain underwater explosive combinations (explosive 
charge weight and standoff), the frequency of this loading can be on the order of the 
natural frequencies of the ship's low frequency flexural modes. Thus, hull flexural 
vibrations, excited by the quasi-periodic bubble pulse loading, can approach a resonant 
response in the low frequency modes. It is clear then that damping can have a significant 
influence on the magnitude of the whipping response of the hull to the bubble pulse 
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Figure 1-2. Effect of damping (with frequency) on the steady state 
vibration of a single degree of freedom system 



2. DAMPING MECHANISMS IMPORTANT IN SHIP WHIPPING 



2.1 Introduction. 

In dynamics, the word damping generally refers to the dissipation of energy during 
the motion of a body or structure This energy eventually is converted to heat by friction 
or generates waves in the surrounding medium When a structure, such as a ship, vibrates 
or moves in on ocean environment, energy can be dissipated within the structure itself, 
through internal components, or externally into the surrounding fluid Specific 
mechanisms of damping, which can be expected to be important in ship and submarine 
whipping, are discussed here and then incorporated into a model for submarine whipping 
in chapter 3. 

2.2 Hull damping. 

The loss of vibrational energy within a complex structure can be accomplished 
through several possible mechanisms, depending on the frequency of the vibration and 
geometry of the structure. In general, low frequency energy (consistent with ship 
whipping) is lost within a structure through friction, which in turn generates heat which is 
conducted, convected and radiated away. This friction can occur among the material 
particles or crystals (within the structural elements) as they move relative to one another 
during material deformation This damping mechanism is referred to as material or 

hysteresis damping. Friction can also occur at the interface between structural elements 

/ 

as they slide relative to one another (i.e. at joints). This damping mechanism is often 
referred to as structural or dry’ friction damping. 

As will be discussed in greater detail subsequently, it is useful to define the hull 
damping in terms of equivalent viscous damping. By doing this, the damping forces 
associated with material and structure can be represented as being proportional to the 
velocity. The resulting equations of motion for the vibration of the hull could then be 
written in the form 

Mx + Cx + Kx = f (2-1) 

where M is the mass matrix of the hull, C is the equivalent viscous damping matrix, K is 
the stiffness matrix, x is the displacement vector, f is the external load vector, and 
superposed dots refer to time derivatives (x is velocity and x is acceleration). Although 
material damping is more representative of the energy dissipation within the material, and 
structural damping more representative of losses in joints, it will be noted that an 
equivalent viscous damping is assumed for most real structural dynamic analyses. The 
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precise nature of the damping is less important than modeling the correct energy loss per 
cycle. 



2.2.1 Material (hysteresis) damping. 

When materials are deformed, energy is absorbed and dissipated by the material 
due to friction during internal reconstruction of the micro and/or macro structure as the 
material deforms. This reconstruction ranges from crystal lattice to molecular scale effects. 
This process, referred to as material damping or hysteresis damping , is discussed in great 
detail by Lazan 3 and Nashif 4 . 

All real materials dissipate energy during cyclic deformation Regardless of the 
precise physical mechanisms involved, energy dissipation can be highly nonlinear, and 
therefore detailed analysis can be very difficult One approach toward quantifying the 
internal damping behavior of real materials is through the hysteresis loop , figure 2-1 . The 
loop plots strain vs. stress during cyclic deformation of a material volume and is obtained 
from experimental measurements. Typically, hysteresis loops representative of 
construction metals such as steel are extremely thin, deviating little from a single line 
(unless the metal is strained well into the plastic range) Thus, for elastic ship hull 
whipping of a steel hull, the damping provided by material hysteresis would be expected to 
be small compared to other damping mechanisms discussed subsequently. 

The area inside the hysteresis loop is the specific damping energy \ D, and is equal 
to the energy dissipated per unit volume of material per cycle 



where a is stress and 8 is strain Equation (2-2) implies that the rate of energy dissipation 
(dD/dt) is proportional to the strain rate. This is equivalent to the hysteresis damping being 
proportional to the strain rate. 

The specific damping energy is very small for most conventional structural 
materials, such as steel, for typical elastic stress-strain levels High-damping alloys and 
viscoelastic materials have higher specific damping energies, and therefore are much better 
at dissipating energy. 



3 Lazan. B. J.. Damping of Materials and Members in Structural Mechanics, Pergamon Press. New York. 
1968 . 

"’Nashif. A.D.. Jones, D.I.G.. and Henderson. J.P.. Vibration Damping. John Wiley and Sons. New York. 
1985 . 




( 2 - 2 ) 
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Figure 2- 1 . Typical hysteresis loops. 

As discussed in the previous section, it is useful to represent hysteresis damping in 
terms of an equivalent viscous damping. An equivalent viscous damping coefficient for 
hysteresis damping can be found by considering the harmonic motion of an equivalent 
spring-viscous damper system (figure 2-2). For this system, the force F(t) necessary to 
cause a displacement x(t) is given by 

F(t) = kx(t) + cx(t) (2-3) 

where k is the spring constant, c is the viscous damping constant, and x is the velocity 
(dx/dt) For harmonic motion of frequency co and amplitude X, the displacement can be 
written in the form 

x(t) = Xsin(cot) (2-4) 

Inserting equation (2-4) into equation (2-3) gives 

F(t) = kXsin(cot) + cXco cos(cot) (2-5a) 

F(t) = kx ±cco -y/x 2 - (Xsincot) 2 (2-5b) 

F(t) = kx ± cco Vx : - x : (2-5c) 

When F(t) is plotted against x(t), equation (2-5c) represents a closed loop similar to a 
hysteresis loop (figure 2-3) The energy dissipated by the damper in one cycle of motion 
is found by substituting equation (2-5a) into equation (2-2): 
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AW = <£ode = J Fdx = J f| — jdt = J (kX sin cot + cXco coscot )(coX coscot )dt =7tcocX : 



dt 



(2-6) 

where AW is the energy loss in the damper in one cycle, which is the equivalent to the 
specific damping energy, D, for hysteretic damping. 




Figure 2-2. Equivalent spring-viscous damper system. 




Figure 2-3. Equivalent spring-viscous damper system. 



For real hysteresis damping, however, it has been found experimentally that the 
specific damping energy, D, is independent of frequency, but approximately proportional 
to the square of the amplitude of the strain 5 . Thus, the equivalent viscous damping 
coefficient, c eq , must be inversely proportional to the frequency, which can be represented 
by 



c 



eq 



_h 

co 



(2-7) 



5 Rao. S.. Mechanical Vibrations. Addison-Wellesley Publishing Co.. Reading. MA. 1990. p. 102. 
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where h is called the hysteresis damping constant As discussed subsequently, the idea 
that the equivalent viscous damping constant for material hysteresis damping is inversely 
proportional to frequency is made use of in the practical modeling of material damping in 
real structures. 

2.2.2 Structural (dry friction) damping. 

In addition to damping within the material of a structure, damping also occurs at 
the interface between structural elements as they slide relative to one another This type 
of damping is usually referred to as structural damping, slip damping. Coulomb damping , 
or simply dry’ friction damping (it amounts to Coulomb or "dry" friction forces). This 
type of damping is considered important when structural elements are joined non-rigidly 
by riveting, clamping, or bolting, where a moderate amount of clamping pressure allows 
for some slip between elements. As for material damping, structural damping is discussed 
in great detail by Lazan 6 , and Nashif 7 . 

Because of the potential geometric complexities, quantitative evaluation of 
structural damping in complex structures is impractical even when numerical techniques 
such as finite element analysis are used However, most structurally significant joints 
making up a ship hull are rigidly welded rather than bolted or riveted. Thus, for elastic 
ship hull whipping, it would be expected that damping provided by structural dry friction 
forces would be very small compared to other damping forces discussed subsequently. It 
should be noted, however, that for certain ship designs (particularly many submarines and 
submersibles), internal decking is often connected to the main structure of the hull through 
various types of sliding deck joints or suspension mechanisms Thus, in some ship 
applications, it is necessary to account for some degree of dry friction damping, 
particularly when considering potentially large amplitude hull flexural motions. 

In general. Coulomb's Law of Dry Friction states that when two bodies are in 
contact, the tangential force required to produce sliding is proportional to the normal 
force acting in the plane of contact. Once sliding has begun, the force remains constant. 
The dry friction force is given by 

f D = UN (2-8) 

where N is the normal force and p is the coefficient of friction which is characteristic for 
the interface of the two surfaces. Thus, the structural dry friction damping force at the 
interface of two sliding surfaces is in a direction opposite to the direction of the relative 



6 Lazan. op. cit. 

7 Nashif. op. at. 
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velocities of the two surfaces but is independent of the displacement or velocity; it 
depends only on the normal force between the two surfaces 

The fact that the dry friction forces remain constant as long as sliding motion exists 
between two surfaces (and are not dependent on displacement or velocity) leads to the 
conclusion that, in each successive cycle (or oscillation) of a structure, the amplitude of 
the motion is reduced by an amount proportional only to the friction forces resulting from 
the sliding within the structure For a simple one dimensional case, the reduction in 
amplitude with each successive oscillation will remain constant as long as relative motion 
exists. 

2.2.3 Modeling material and structural damping in complex structures. 

If the material and structural damping properties of a structure could be 
quantitatively determined throughout the structure, a finite element procedure could be 
used to include the distributed damping properties (as element damping matrices) in the 
complete dynamic solution of the structural motions. In practice, however, direct 
evaluation of distributed structural damping properties in complex structures does not 
lend itself well to the numerical solutions for structural dynamics (although distributed 
damping has been included in some simple finite element applications 8 ). Material and 
structural damping (as previously defined) are typically expressed in terms of damping 
ratios established from experiments, rather than by means of an explicit damping matrix. In 
fact, for many applications, there is no need to express the damping explicitly by means of 
a damping matrix, but rather only in terms of the modal damping ratios (£ n ) which can be 
applied through the principle of modal superposition. This concept is often called modal 
damping. An extension of modal damping in which an explicit damping matrix is defined 
through a linear combination of the mass and stiffness matrices is called Rayleigh 
damping. 

Modal damping and modal superposition. The method of dynamic analysis 
using modal superposition is discussed in numerous texts (Clough and Penzien 9 , Smith 10 , 
for example) and is discussed briefly here only for continuity. Fundamentally, the original 
set of N coupled equations of motion (for an N degree of freedom system) are 
transformed to a set of N independent normal (modal) coordinate equations. The 

8 Kaliske. M.. Jagusch. J.. Gebbeken. N . and Rothert. H . "Damping models in finite element 
computations". Structural Dynamics-Proceedings of the 2nd European Conference on Structural 
Dynamics (I 'ol. I). Rotterdam. 1993, pp. 585-591. 

9 Clough. R. W. and Penzien, J.. Dynamics of Structures (second edition), McGraw-Hill. New York. 1993. 

10 Smith. J.W.. I ’ ihration of Structures - Applications in Civil Engineering Design. Chapman and Hall. 
London. 1988. 
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dynamic response can be obtained by solving separately for the response of each normal 
(modal) coordinate and then superposing these to obtain the response in the original 
coordinate system. 

The equations of motion for a multiple degree of freedom dynamic systems can be 
written in matrix form: 



Mx + Cx + Kx = f (2-9) 

where, in general, matrices M, C, and K all have non-diagonal terms (i.e. the equations of 
motion are coupled) The modal equations of motion are obtained from equations (2-9) 
by applying a modal transformation, 

x = <t>q (2-10) 

where O is the matrix whose columns are the mode shape vectors 

0 = [(t),(t) ; ...<j) N ] (2-11) 

(<j)j is the vector representing the shape of the ith mode), and q is a vector of normal or 
modal coordinates. Pre-multiplying each term by the transform of the /rth mode shape 
vector, yields an equation of motion for each mode (uncoupled from other modes) of the 
form: 

MA+C„q n +K„q„=Q„ (2-12) 

where: 

M n = 4> n r M<h = modal mass coefficient 

C n = <}) n r C<j) n = modal viscous damping coefficient 

K n = <l>nK<t> n smodal stiffness coefficient 

Q n = <j> n r f s modal force (2-13) 

Equations (2-13) for modal mass and modal stiffness come about because of the 
orthogonality property of the normal mode shapes, which can be extended to modal 
viscous damping if the orthogonality condition is also assumed to apply to the damping 
matrix (as will be discussed subsequently). If equation (2-12) is divided by the modal 
mass coefficient, the modal equations of motion may be written in an alternate form: 



q'n +2Cn® n q n +®nq n = 



On. 

M. 



(2-14) 



where to n is the (modal) natural frequency, and £ n is defined as the modal viscous 
damping ratio 



c„ 



c n 

2o n M n 



(2-15) 
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The (modal) natural frequencies (co n ) and mode shape vectors (<J> n ) are found by solving 
the generalized eigenproblem for the undamped free vibration of the dynamic system: 

K<l> n =o^M4> n (2-16) 

As discussed previously, it is often more convenient (and physically reasonable) to 
define the damping of a complex multiple degree of freedom system using a viscous 
damping ratio for each mode, rather than to try to explicitly evaluate coefficients of the 
damping matrix C, because the modal damping ratios can be determined experimentally 
(or estimated with adequate precision in many cases). 

With the (uncoupled) modal equations of motion defined, the modal displacements 
can be found, and the total displacement can be found by summing the modal 
contributions using equation (2-10): 

x = Z4>„q n (2-17) 

n=l 

Thus, the total response of a system modeled using modal superposition can be thought of 
as a sum of the responses of N independent single degree of freedom systems. 

Rayleigh (proportional) damping. Under some circumstances, it is desirable or 
practicable to use modal superposition to uncouple the equations of motion in solving for 
the dynamic response of a system. In other circumstances, it may be desirable to use the 
basic concept of modal superposition, but develop an explicit damping matrix (C) which 
may be used to solve the complete equations of motion in the time domain. This can be 
done by applying what has become known as Rayleigh damping or proportional damping 
(named after Lord Rayleigh, who first suggested its use). Rayleigh or proportional 
damping is discussed in varying detail in numerous sources (Clough and Penzien 11 , 

Bathe 12 , for example). 

Under Rayleigh or proportional damping, it is assumed that the damping matrix 
can be adequately represented by making it a linear superposition of the mass matrix (mass 
proportional damping), the stiffness matrix (stiffness proportional damping), or both. By 
applying a linear superposition, the orthogonality condition which applies to the mass and 
stiffness matrices (as discussed previously) would also apply to the damping matrix. Thus, 
the damping matrix could be written in the form: 

C = aM + pK (2-18) 



"Clough and Penzien, op. cit. 

l2 Bathe, K.J., Finite Element Procedures in Engineering Analysis, Prentice-Hall, Englewood Cliffs, NJ, 
1982 . 
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where a and (3 are proportionality constants The damping associated with equation (2-18) 
may be recognized by evaluating the generalized modal damping associated with it by 
combining it with equations (2-13) and (2-15), and solving for the modal damping ratio 

C n = <l>nC<t> n = ct>n[aM +pK]<t> n = aM n + PK n (2-19) 



Cn = 



2co 

n n 



a ( Pco n 



2o, 



( 2 - 20 ) 



Thus, the damping associated with any particular mode can be thought of as a linear 
combination of a factor which is inversely proportional to frequency (the mass 
proportional term), and a factor which is directly proportional to frequency (the stiffness 
proportional term). An important point to note is that the Rayleigh damping model is 
consistent with the discussion of material (hysteresis) and structural damping presented 
earlier. Thus, Rayleigh damping can be thought of as providing a good model to account 
for material and structural damping (for elastic whipping, as discussed in the previous 
sections). 

The relationship between damping ratio and frequency for Rayleigh damping can 
be illustrated as shown in figure 2-4. The cases of purely mass proportional damping 
(P = 0) and purely stiffness proportional damping (a = 0) are shown separately from the 
combined or general case. 

As illustrated in figure 2-4, the two Rayleigh damping proportionality constants, a 
and P, can be evaluated by solving two simultaneous equations if the (modal) damping 
ratios, C a and £ b , associated with two specific (modal) frequencies, co a and (o b , are 
known (i e. from experiment). Substitution of the two damping ratios and frequencies into 
equation (2-20), combining the resulting equations into matrix form, and carrying out a 
matrix inversion gives the solution of the two damping proportionality constants: 




1 /( 0 . 

l/co, 



(0 



oo. 



( 2 - 21 ) 



= 2 - 
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- 1 /( 0 , 



1 / (0 




(2-22) 



It is apparent from figure 2-4 that the use of Rayleigh damping provides for very 
high damping ratios for modes of very high frequencies (<a much greater than o b ). The 
end result is that the responses of very high frequency modes are effectively eliminated by 
the high damping imposed by the Rayleigh damping. In the case of elastic ship whipping, 
where only the lower few modes are of concern (as discussed previously), this limitation 
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can be considered acceptable as long as the higher frequency used for evaluation of the 
Rayleigh proportionality constants (to b ) is set among the higher modes expected to 
contribute significantly to the whipping response 



Damping ratio 




Figure 2-4. Relationship between damping ratio and frequency for Rayleigh damping. 

2.3 Internal damping. 

Within any complex structure, such as a ship, energy can be lost (or simply 
transferred) to "separate" dynamic systems which are in some form "attached" to the main 
structure. For ships, two broad categories could be considered under this category. 
Discrete ititemcil masses or structures such as machinery and equipment connected 
indirectly to the hull through decking, machinery foundations, and shock or sound mounts, 
can act as dynamic subsystems, redistributing and absorbing dynamic energy transferred 
from the hull. Additionally, onboard liquids which have free surface can act as dynamic 
subsystems, redistributing and absorbing dynamic energy by sloshing. Because of the 
large number of discrete internal structures and liquid storage tanks found onboard most 
ships, it is important to consider the mechanisms through which energy can be transferred 
and absorbed into these internal damping devices 

2.3.1 Damping through discrete internal structures. 

Heavy machinery, electronic equipment, and other discrete masses onboard a ship 
can be considered to have a significant potential for dissipating hull whipping energy. 
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Figure 2-5 illustrates notional discrete masses, connected to the hull via mounts, internal 
decking, foundations, etc. As will be discussed subsequently, the dynamic response of 
more complex internal structures could be modeled as easily as individual masses by 
applying the principle of modal superposition (under the assumption of linear response). 
For purposes of this analysis, the term "discrete internal structure" will be used to refer to 
any mass or discrete structure that may be modeled as a discrete internal dynamic system 
in terms of transfer and absorption of hull whipping energy. While it is certainly not 
desired that heavy machinery or electronic equipment absorb large amounts of energy, it 
must be assumed that their ultimate connection to the hull provides the potential for the 
transfer of energy from the whipping hull. 




Figure 2-5. Discrete masses connected to the hull via 
mounts, internal decking, foundations, etc. 

Whether or not discrete internal structures connected to a hull provide damping to 
hull whipping depends upon the sizes of the masses, the dynamic properties of the 
mechanisms through which they are connected to the hull (decking, foundations, mounts, 
etc.), the dynamic properties of the hull, and the frequency of hull whipping (frequency of 
excitation). In section 1.2, it was stated that the controlling mechanism in the response of 
a system to excitation (stiffness, damping, or mass) is determined by the relation of the 
excitation frequency (or frequencies) to the natural frequency (or frequencies) of the 
system. The maximum amount of energy is transferred when the frequency of excitation is 
equal to a natural frequency (at resonance) In the case of discrete internal structures 
connected to a whipping hull, it is evident that maximum energy would be transferred 
when the frequency of hull whipping is equal to a natural frequency of the system 
composed of the hull, the discrete internal structure, and the connection mechanisms. An 
understanding of the dynamic effects of discrete internal structures connected to a 
whipping ship can be made by first considering the effects of a single mass connected to a 
(rigid) vibrating structure, where the structure and mass are free to vibrate in only one 
direction (a two degree of freedom system), and then, using the principle of modal 
superposition, extending this to more complex arrangements of multiple masses connected 
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to a vibrating structure, such as a whipping ship hull (i.e. multiple degree of freedom 
systems). 

The two degree of freedom system (dynamic vibration absorber) A single 
mass connected to a much larger vibrating structure can be used to dissipate or absorb 
vibrational energy of the larger structure. In machinery and structural dynamics, this 
concept has gained practical significance in the application of the dynamic vibration 
absorber The concept of the dynamic vibration absorber is presented by Den Hartog 13 , 
and their design and application is discussed in detail in numerous other sources 
(Korenev 14 , Hunt 15 , for example). A small auxiliary mass is connected to a harmonically 
vibrating structure by a spring of specified stiffness (the location of the mass is selected to 
have the maximum effect on the significant resonant mode of the vibrating structure) In 
the case where damping is provided by the connecting mechanism (e.g. by a dashpot), the 
auxiliary mass is called a damped dynamic vibration absorber or tuned mass damper 
Treating the main structure and vibration absorber as a two degree of freedom dynamic 
system (figure 2-6), Den Hartog showed that the amplitude of harmonic vibration of the 
main structure could be reduced to zero at a specified excitation frequency (or reduced 
over a range of frequencies). This could be accomplished by "tuning" the properties of 
the vibration absorber (adjusting the mass, stiffness of the spring, and damping in the 
dashpot). The downside in the application of dynamic vibration absorbers appears in the 
potentially large amplitude response (displacement and acceleration) of the absorber 
masses. In the case of heavy machinery and electronic equipment onboard a whipping 
ship, this could translate to high displacements and accelerations on the delicate 
equipment, and large forces on the mounts and decking. 

For the simple dynamic vibration absorber attached to a large vibrating structure, 
the equations of motion for this simple two degree of freedom system can be derived by 
considering dynamic equilibrium of each of the masses separately. Figure 2-7 illustrates 
each of the masses and the dynamic forces acting upon them. From equilibrium, the 
equations of motion may be written: 

MX + KX = F ext +c(x - X) + k(x - X) (for the main mass) (2-23) 

mx +c(x - X) + k(x - X) = 0 (for the absorber mass) (2-24) 

where capital letters refer to the main structure and lower case letters refer to the 
absorber, as shown in figure 2-6 These equations of motion could be written in an 
alternate form by defining the force exerted on the main structure by the absorber : 

13 Den Hartog, J P., Mechanical I'ibrations. (4th edition), McGraw-Hill. New York. 1956. 
l4 Korenev and Re/.nikov. Dynamic l ibration Absorbers . John Wiley and Sons. New York. 1993. 

15 Hunt. J.B.. Dynamic Vibration Absorbers. Mechanical Engineering Publications. London. 1979. 
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Berber = k( X - X) + C( X - X) (2-25) 

which, of course is the negative of the force exerted on the absorber by the main structure. 
Thus, the equations of motion could be written: 

MX + CX + KX = F ext + F abs (for the main mass) (2-26) 

Fabsorber = c(x - X) + k(x - X) =r-mx (for the absorber mass) (2-27) 
Equations (2-26) and (2-27) form a set of 2 simultaneous differential equations (coupled) 
which may be solved under various conditions for displacements, velocities and 
accelerations. 



I 1 




Main 

Structure 



Figure 2-6. Dynamic vibration absorber. 




1 = 



X<t) 




k<x-X> c(x-X) 



ABSORBER MASS 
Figure 2-7. Dynamic forces on each mass of a dynamic vibration absorber. 



Multiple degree of freedom systems. For more complex arrangements of 
internal structure, where many degrees of freedom would be required to properly model 
the dynamic response, the displacements, accelerations, and velocities would be vectors, 
and a larger number of simultaneous differential equations would be required. For 
multiple degree of freedom main and internal structures, the equations of motion could be 
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written just as equations (2-23) and (2-24), but in matrix form. Thus, there would be one 
simultaneous equation for each degree of freedom. 

An alternative approach for multiple degree of freedom systems (particularly 
where the internal structure is complex and a large number of degrees of freedom would 
be required to adequately represent the dynamics), would be to use the principle of modal 
superposition under the assumption of linear response (as discussed in section 2 2.3) to 
approximate the response of the internal structure By including only those resonant 
response modes of the main structure and the internal structure which would contribute 
significantly to the overall response, the complexity of approximating the response of the 
main structure is greatly reduced. Figure 2-8 depicts a notional complex internal 
arrangement of discrete masses with many degrees of freedom. 



Using modal superposition, an equivalent mechanical model of the dynamic 
behavior of the discrete internal structure may be made as shown in figure 2-9. Each 
mode of resonant response of the discrete internal structure is represented by an 
equivalent modal mass, modal damping and modal stiffness (mass m 0 represents that 
portion of the discrete internal structure responding as a rigid body with the main 
structure). A variation of the dynamic equilibrium equation (2-24) for the dynamic 
vibration absorber can be written using the principle of modal superposition. The equation 
of motion for the /rth mode of the internal structure can be written: 



where x n is the /?th modal displacement of the internal structure (equivalent to q n in 
section 2.2.3) and X is the displacement of the main structure (hull). 

The equations of motion for a system made up of an internal structure attached to 
a large vibrating structure (the hull) can be derived by considering dynamic equilibrium of 




Figure 2-8. Complex internal arrangement of discrete masses with 
many degrees of freedom 




(2-28) 
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the main structure and the superposed modal masses of the internal structure. The forces 
applied to the main structure by the internal structure can be obtained by modal 
superposition. Figure 2-10 illustrates the dynamic equilibrium for the main structure. 




structure 



Figure 2-9. Equivalent mechanical model of the dynamic behavior of 
a discrete internal structure (using modal superposition). 




Figure 2-10. Dynamic equilibrium for the main structure. 



From dynamic equilibrium, the equation ot motion for the main structure can be 

written 

[MX + CX + KX] = F ext + F mt (2-29) 
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where F ext is an external force applied to the main structure, and F mt is defined as the force 
exerted on the main structure by the internal structure. F int can be determined by 
considering dynamic equilibrium of the internal structure using modal superposition, as 
shown in figure 2-11. By dynamic equilibrium of the internal structure: 

F ,nt =-i> n X n =X[ C nU„ -X) + k n (x n -X)l (2-30) 

n=0 n~0 



E mx 



jx(t) 



T E k(x-x) T E c(s-x) 

Figure 2-11. Dynamic equilibrium of the internal structure using modal superposition. 



It should be noted that the internal structure rigid body case (n = 0) can be taken out of 
equation (2-30) by noting: 



E[ c n(’ i „-X) + k„(x n -X)] = -^m n x n =-m 0 X-£ 



m x = -m, 



n=0 



|X + E[ C n(^n-X) + k n (.X n -X)] 



(2-31) 

With this simplification, the force exerted on the main structure by the internal structure 
can be written: 



F,=-i».Xtlk(j,-X)tk.(x,-X)l (2-32) 

n=l 

Equations (2-29), (2-28), and (2-32) form a set ofA/W simultaneous equations (coupled) 
which can be solved under various conditions for displacements, velocities and 
accelerations 

It should be noted that, for cases where more than one degree of freedom is 
important for the main structure (e.g. the main structure moves horizontally as well as 
vertically), the above method could be extended considering the modal response of the 
internal structure in each of the degrees of freedom of the main structure. Thus, the main 
structure displacement would be represented by a vector (X), the internal structure modal 
displacement would be represented by a vector (x n ), and the forces (F mt and F ext ) would 
be represented by vectors. 
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An application of damping through discrete internal structures will be included in 
the development of a model for submarine whipping in chapter 3. 

2.3.2 Damping through sloshing liquids. 

Liquid sloshing has been the focus of considerable attention in the space industry, 
due primarily to concerns about fuel-rocket interaction forces, and recognition that liquid 
sloshing resonance can result in significant dissipation of energy and can, therefore, be 
capitalized upon to suppress vibrations in large structures 16 . Similarly, methods to 
account for the effects of sloshing liquids on the natural frequencies and damping of 
offshore platforms have been developed within the marine industry 17 The occurrence of 
resonance, effects of nonlinearities, and contribution of devices such as baffles have been 
the object of numerous other studies. It can be projected that liquids in tanks onboard a 
whipping ship could undergo sloshing resonance, resulting in the damping of whipping 
energy. Fundamentally, any onboard liquid with free surface has the potential for 
absorbing whipping energy through liquid sloshing (liquid without free surface acts only as 
a rigid mass). Figure 2-12 illustrates a notional liquid storage tank, restrained to move 
with the whipping hull, and the resulting sloshing (or standing waves) generated due to the 
existence of the free surface of the liquid in the tank. 



The liquid in a moving tank can be expected to respond in a variety of ways 
dependent upon many variables (tank geometry, directions of tank motion, amplitudes and 
frequencies of tank motion, properties of the liquid, etc ). Translational or pitching 



16 Welt. F., and Modi. V.J.. "Vibration damping through liquid sloshing". Diagnostics, Vehicle Dynamics 
and Special Topics: 1 2th Biennial Conference on Mechanical I 'ihrations and Noise. ASME. 1989. 
17 Vandiver. J.K.. and Mitome. S.. "Effect of liquid storage tanks on the dynamic response of offshore 
platforms". Dynamic Analysis of Offshore Structures: Recent Developments. CML Publications. 
Southampton. 1982. 




x(t) 



Figure 2-12. Tank with sloshing liquid. 
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motions of a tank lead primarily to lateral sloshing (usually antisymmetric liquid mode 
shapes). Tank motions normal to the equilibrium free surface lead primarily to vertical 
sloshing (usually symmetric liquid mode shapes). Under conditions of abrupt or sudden 
acceleration, liquid impact with tank overheads and opposing bulkheads may provide 
significant loading on tank boundaries. The high accelerations associated with a whipping 
ship may require consideration of any or all of these to account for the transfer of 
whipping energy. 

Analysis of liquid sloshing can be carried out by numerous means. Researchers in 
the field of liquid sloshing have used several principle methods to determine the important 
dynamic characteristics desired for various applications 18 . In many applications, potential 
flow models, using linearized and nonlinear free surface conditions, have been developed 
and used to approximate the fluid motion in various container shapes (Bauer 19 , Hutton 20 , 
Welt and Modi 21 , for example). In order to properly account for energy transfer and/or 
dissipation, incorporation of liquid impact loading and viscous effects has been 
accomplished using equivalent mechanical models in conjunction with potential flow 
models (Dalzell 22 , Vandiver and Mitome 23 , Bauer 24 , for example). 

Potential flow model. Fundamentally, lateral and vertical sloshing can be thought 
of as the motion of standing waves within a moving rigid container or tank. Under the 
assumptions of incompressible and irrotational flow, the motion of the liquid in a tank can 
be approximated well by solving the potential flow problem for the wave motion subject to 
appropriate boundary conditions (all variables are functions of time). The velocity 
potential function <t> (defined by v = V O where v is the fluid velocity vector and V is the 
vector differential operator) represents solution of the differential equation (Laplace 
equation) 

V : <t>-0 (within the fluid domain) (2-33) 

subject to the following boundary conditions: 



18 Abramson, H.N. (ed.), The Dynamic Behavior of Liquids m Moving Containers, .VASA SP-106. 1966. 
19 Bauer. H F.. "Theory of the Fluid Oscillations in a Circular Ring Tank Partially Filled with Liquid". 
.VASA Technical Note D-557. 1960. 

20 Hutton. R E., "An Investigation of Resonant. Nonlinear. Non-Planar Free Surface Oscillations of a 
Fluid". NASA Technical .Vote D- 1 870. 1963. 

21 Welt and Modi. op. cit. 

22 Dalzell. J.F.. "Liquid Impact on Tank Bulkheads". The Dynamic Behavior of Liquids in Moving 
Containers, XAS. I SP-106. H.N. Abramson (ed.). 1966. pp. 353-372. 

23 Vandiver and Mitome. op. cit. 

24 Bauer, H F.. "Nonlinear Mechanical Model for the Description of Propellant Sloshing". Al.-l i Journal. 
Vol. 4. No. 9. 1966. pp. 1662-1668. 



29 



(2-34a) 



rO 

= V n (on the wall of the tank) 

cn 



c<b rO f<J> 1 

— - + g + 2 VO- V + -VO- V(VO- VO) = 0 

ct" cy ct 2 

(at y = r| (on the free surface)) (2-34b) 



where n is the normal vector of the tank wall, V n is the velocity of the tank wall normal to 
its surface, g is the acceleration due to gravity, and q is the free surface elevation The 
free surface boundary condition, equation (2-34b), represents the complete nonlinear free 
surface boundary condition 25 . The linearized free surface boundary condition (linearized 
about the mean free surface) would be written: 



c~ O cO . . 

-z^- + g— = 0 (at y = 0) (2-34c) 

ct' cy 

The fluid velocity vector (v) can be obtained from (by definition of the velocity potential): 

v = Vd> (2-35) 

The hydrodynamic pressure (P h ) exerted on the tank wall by the sloshing liquid can be 
obtained from the unsteady Bernoulli equations: 



-P 



cO 

ct 



+ -|V<D 
2 



(2-36) 



where p is the liquid density. Finally, the hydrodynamic force exerted on the tank wall by 
the sloshing liquid can be obtained by integration of the hydrodynamic pressure over the 
area of the tank wall: 



F h =jj(P h -n)dS (2-37) 

S 

where F is the hydrodynamic force vector (on the tank wall) and S is the surface of the 
tank wall (below the free surface). Thus, for a specified tank geometry and motion, the 
hydrodynamic forces exerted by the sloshing liquid on the tank can be approximated using 
a potential flow model. 

Equivalent mechanical model. As mentioned previously, in order to properly 
account for energy dissipation, viscous effects at the tank wall must be included. While it 
is difficult to include viscous effects in the potential flow solution, viscous damping may 
be included in an equivalent mechanical model for the sloshing. In an equivalent 

25 Ne\vman. J.N.. Marine Hydrodynamics. The MIT Press. Cambridge. MA. 1986. p. 247. 
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mechanical model, each of the sloshing modes (or modes of standing wave motion) is 
represented by an equivalent spring-mass-damper as shown in figure 2-13. In typical 
practice, linear viscous damping is usually assumed for the (modal) dampers in the 
mechanical model- 6 , however, nonlinear mechanical modeling techniques do exist 27 



X(t) 




Figure 2-13. Equivalent mechanical model of a tank with sloshing liquid 
using modal superposition. 



In an equivalent mechanical model, the size of the equivalent masses for each of 
the resonant sloshing modes (m, , m, ,...) may be determined by requiring that the 
(hydrodynamic) force exerted on the tank wall by the equivalent mass-spring system, for 
the case of zero damping, be equal to the force calculated from the potential flow theory 
for each resonant sloshing mode (m 0 is the portion of the liquid in the tank which acts as a 
rigid body). 



f n = m n* n => m„ = -^~ 

X. 



(2-38) 



- Silverman. S.. and Abramson. H.N.. "Damping of Liquid Motions and Lateral Sloshing". The Dynamic 
Behavior of Liquids in Moving Containers, \ T ASA SP-106. H.N. Abramson (ed.). 1966. pp. 105-144. 
27 Bauer. op. cit. 
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where f n h is the hydrodynamic force exerted on the tank wall by the mh sloshing mode 
The equivalent modal stiffnesses (k, , k, ,...) may be determined from the natural 
frequencies (co n ) calculated from the potential flow theory for each resonant sloshing 
mode: 



© 




© :m„ 



(2-39) 



Damping coefficients for each resonant sloshing mode can be determined 
experimentally for a given tank configuration via a logarithmic decrement method (section 
1.2). The modal damping ratio for the mh mode (£ n ) can be determined from the log 
decrement (5) measured for each mode by 



C, n = , . .. = w — (for small damping) (2-40) 

V(27r) 2 +5 3 2 tc 

where the log decrement is defined by equation (1-1). The modal damping ratio is the 
ratio of the modal damping coefficient to the critical modal damping coefficient (c/ c .), 
where the critical modal damping coefficient represents that value of damping where the 
motion would not oscillate, but would be non-periodic. The modal damping coefficients 
( c, , c ; , . . . ) can be determined from 

c„ =2<;„m n <0, (2-11) 

Referring to figure 2-13, the equation of motion for the mh equivalent mass can be 
written in terms of the displacement, velocity, and the acceleration of the equivalent mass 
as follows: 

m,S n +c n (x n -X)+k,(x,-X) = 0 (2-42) 

where X is the displacement of the tank, x n is the displacement of the mh sloshing mode 
(equivalent mass), and velocities and accelerations are represented by superposed dots. 
The total hydrodynamic force on the wall of the tank can be found by summing the forces 
contributed by each mode: 



Fh =-Z m n*n =- m oX-X m n*„ (2-43) 

n~0 n-1 

Similar to the discrete internal structure, each of the resonant modes of liquid sloshing 
can be applied to the main structure (the hull). 

As an example, consider the lateral motion of a rigid structure with an integral 
laterally sloshing tank. The equivalent mechanical model could be as shown in figure 
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2-14 The equation of motion for each of the equivalent liquid masses (modal masses) is 
given by equation (2-42), and the total hydrodynamic force on the wall of the tank due to 
the liquid sloshing is given by equation (2-43) The equation of motion for the main 
structure can then be found by adding the total hydrodynamic force into the equation of 
motion for the structure without the liquid: 

[ MX + C X + KX ] = F ext + F h (2-44) 

where F ext is an external force applied to the main structure. The hydrodynamic force (F h ) 
may be written (following the same principles of modal superposition used for the internal 
structure): 

=-2>. s . * -m„X+£[c.<*. -X) + M*. -X)] (2-45) 

n-0 n-1 

It can also be noted that m 0 , the portion of the liquid which acts as a rigid body, can be 
found from. 

M i = Z m n => m o = M, -f]m n (2-46) 

n-0 n=l 

where M, is the total mass of liquid in the tank. 




Figure 2-14. Equivalent mechanical model for the lateral motion of a rigid 
structure with an integral laterally sloshing tank. 



In terms of real applications, such as in the whipping of a ship hull, equation (2-45) 
may be truncated because the equivalent masses m n and modal displacements x a grow 
smaller with n (the rate at which they grow smaller depends upon the geometry of the 
tank). Ali terms should be included for which the natural frequencies of the tank (® n ) are 
less than or equal to the natural flexural frequency of the overall structure (less than the 
fundamental flexural/whipping natural frequency). If equation (2-45) is truncated to 
include A'’ modes, then equations (2-44), (2-45), and (2-46) form a set of N- 1 
simultaneous equations of motion for the total dynamic system consisting of the main 
structure and the sloshing liquid. 

It should be noted that, for cases where more than one degree of freedom is 
important for the main structure (e g. the main structure moves horizontally as well as 
vertically), the above method could be extended considering the modal response of the 
sloshing liquid in each of the degrees of freedom of the main structure Thus, the main 
structure displacement would be represented by a vector (X), the sloshing liquid modal 
displacement would be represented by a vector (x n ), and the forces (F h and F exl ) would be 
represented by vectors. 

An application of liquid sloshing will be included in the development of a model for 
submarine whipping in chapter 3. 

2.4 External (fluid) damping. 

The loss of vibrational energy of a whipping ship into the surrounding fluid 
medium can be accomplished through several general mechanisms. If the ship is 
sufficiently close to the free surface, or its motion is of sufficient velocity and/or 
frequency, surface and pressure waves of non-negligible energy can be generated. This 
form of fluid damping is often referred to as wave radiation damping. Additionally, 
whipping energy can be lost through mechanisms associated with the viscosity of the 
surrounding fluid. These mechanisms can be lumped under the general heading of viscous 
fluid damping , but include specific mechanisms associated with time-dependent viscous 
boundary layer shear forces ("skin friction" drag) as well as boundary layer separation 
("form" drag). 

In general, the dominant hydrodynamic forces on a structure are determined 
primarily by the characteristics of the fluid flow around the structure For a whipping ship, 
the relative fluid flow around the hull would be generally oscillatory in nature. For 
oscillatory flow, the following non-dimensional parameters are often used to characterize 
the flow: 



34 



Peak Reynolds' number, Rn = 



U 0 D 



Keulegan-Carpenter number, KC = 



U 0 T 

D 



Viscous-frequencv parameter, B = 

KC vT 

where U is the characteristic relative free stream velocity (U 0 is the amplitude of the 
characteristic relative free stream velocity or maximum relative velocity), D is the 
characteristic dimension of the structure (i.e. the diameter), v is the kinematic viscosity of 
the fluid, and T is the characteristic period of oscillation of the relative fluid flow. 
Additionally, the flow can be influenced significantly by the geometry of the structure, free 
surface effects, sea floor effects, roughness on the surface of the structure, and the 
orientation of the structure relative to the flow It may also be shown that for 
harmonically oscillating flow around a fixed circular cylinder, KC = 2xA / D, where A is 
the amplitude of the oscillation of the fluid particle displacement. Thus, KC may be 
thought of as representing an amplitude of motion betw een a structure and the fluid 
relative the characteristic dimension of the body. 

These parameters, in addition to the other influencing factors listed above, 
determine the characteristics of the flow, and thus the major components of the 
hydrodynamic force on the structure. As will be discussed in greater detail subsequently, 
these characteristic parameters (i.e. the characteristics of the flow) effect the occurrence of 
boundary layer separation, and thus a major component of viscous fluid forces (as will be 
discussed subsequently, boundary layer separation forces tend to dominate over boundary 
layer shear forces for bluff bodies). For fully submerged bodies, if KC is large (i.e. U 0 is 
large and D is small), then boundary layer separation is likely to occur, and the major 
hydrodynamic forces are likely to be associated with viscous drag effects (the so-called 
"drag-dominated regime"). If both Rn and KC are large (i.e. U 0 is large but D is not 
small), then the added forces associated with the inertia of the surrounding fluid (called 
fluid inertia forces) must also be considered in addition to the drag forces (the so-called 
"Morison regime"). If KC is small (Rn may still be moderately large), then boundary 
layer separation is not likely to occur and the only viscous forces would be due to 
boundary layer shear. If, in addition to KC being small, the structure is on or near a free 
surface then radiation of surface waves tends to dominate the hydrodynamic damping (the 
so-called "diffraction/radiation" regime). It is important to note that, for complex 
structures such as ships, different portions of the structure may be subject to different 
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types of hydrodynamic loading (depending upon the local characteristics of the structure 
and the relative fluid flow) 



2.4.1 Wave radiation damping. 

Fundamentally, waves are radiated whenever a structure moves in a fluid medium 
As mentioned above, when KC is small (i.e. when D is large), and the structure is on or 
near a free surface, then forces associated with wave radiation are generally important. In 
typical practice, panel or boundary element methods are used to solve the boundary-value 
problem in analyzing hydrodynamic forces due to incident, diffracted and radiated 
pressures. Boundary element methods are based upon potential flow theory, which 
neglects the effects of viscosity. However, there has been some effort in recent years to 
include some of the effects of viscosity within the boundary element method (as will be 
discussed subsequently) 

Potential flow theory and its application to the hydrodynamics of structures is 
discussed in great detail in numerous sources (Newman 28 , Chakrabarti 29 , Faltinsen 30 , for 
example). Under potential flow theory (assuming inviscid, irrotational flow), the fluid 
velocity within the fluid domain may be represented as the gradient of a scalar potential 
function 



where v is the fluid velocity vector, 0(x,y,z,t) is the scalar velocity potential function 
(written here in terms of Cartesian coordinates), and V is the vector differential operator. 
The velocity potential must satisfy Laplace's equation 



within the fluid domain, and satisfy boundary conditions specific to the geometry of the 
application. For the case of a body of general shape moving near a free surface in a fluid 
of finite depth (figure 2-15), the velocity potential must satisfy the following boundary 
conditions: 



28 Nc\vman. op. cit. 

29 Chakrabarti. S.K.. Hydrodynamics of Offshore Structures. Springcr-Verlag. London. 1987. 
30 Faltinscn. O.M.. Sea Loads on Ships and Offshore Structures. Cambridge University Press. UK. 1990. 




(2-47) 




^ 2 -> 2 a 2 

ox cy cz 



(2-48) 



= V • n (on the surface of the body) 
on 



(2-49a) 
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(2-49b) 



- — = 0 (at y = -h (for finite water depth)) 

cy 




(at y = r| (on the free surface)) 



(2-49c) 



where n is the normal vector of the body surface (pointing from the body into the fluid 
domain), V is the velocity vector of the body, g is the acceleration due to gravity, r| is the 
free surface elevation above the mean, and h is the mean depth of the water. The free 
surface boundary condition, equation (2-49c), represents the complete nonlinear free 
surface boundary condition 31 . The linearized free surface boundary condition (linearized 
about the mean free surface) would be written: 




(2-49d) 




h 




7 7 7 7 7 7 7 7 7 7 7 7 7 



Figure 2-15. Body of general shape moving in fluid with free surface. 



3 'Newman, op. cit.. p. 247. 
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Additionally, a radiation condition must be established (as a boundary condition limit at 
infinity) in order to ensure that radiated waves propagate away from the body. For 
solution of equation (2-48) in the time domain, initial conditions are also required' 

<t> -> 0 



c<t> 

» 0 (as t — > -oo ) 

c t 



(2-50) 



Solution of Laplace’s equation (2-48) for the velocity potential, <t>, is obtained by the 
method of integral equations using Green's theorem which can be written in the general 
form: 



| (<t>V : G - GV : <t>)dV = J (cfiVG - GV<t>)- ndS (2-51) 

V S 

where V is defined as the volume bound by surface S, where S includes the body surface, 
the free surface, the bottom surface, and the surface at "infinity" (figure 2-16) An 
appropriate Green function (G) can be obtained which satisfactorily represents the source 
potentials and their derivatives for various boundary condition combinations (Wehausen 
and Laitone 32 , Newman 33 , for example). 

Once the velocity potential function is obtained, the fluid velocity vector in the 
fluid domain can be found by definition of the potential function as written in equation 
(2-47). The hydrodynamic pressure exerted on the body by the fluid is obtained from the 
unsteady Bernoulli equations: 



P = -p 



c<t> 1 | 

+ — VO 

a 2 



(2-52) 



where P is the hydrodynamic pressure (a function of position) and p is the fluid density. 
The hydrodynamic force exerted on the body (or the panel of interest) by the fluid is 
obtained by integrating the hydrodynamic pressure over the submerged portion of the 
surface: 



F = JJ(P-n)dS (2-53) 

S 

where F is the force vector on the surface or panel, S is the submerged surface of the body 
or panel, and n is the direction normal of the panel surface. 



32 Wehausen. J.V., and Laitone. E.V.. "Surface Waves”. Encyclopedia of Physics (Vol. 9). Springer- 
Verlag. Berlin. 1960. 

33 Ne\vman. J.N.. "Algorithms for the free surface Green function". Journal of Engineering Mathematics. 
19, 1985. pp. 57-67. 
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Figure 2-16. Surfaces of integration for Green's theorem. 



For slender bodies (where length is much greater than width), application of 
two-dimensional diffraction/radiation potential theory has proven accurate and useful, 
particularly at higher frequencies 34 . Several two-dimensional approaches are available in 
the literature, the most common of which is strip theory ,35 . In strip theory, the slender 
body is divided into many thin slices or strips, each of which has approximately constant 
cross-section, so that two-dimensional potential theory could be applied to determine the 
approximate hydrodynamic forces on each strip independently (figure 2-17). For rigid 
body hydrodynamics applications, the hydrodynamic forces on each strip could then be 
integrated along the length of the body to determine the total hydrodynamic force on the 
body. Alternatively, similar to panel methods, the hydrodynamic forces on each "strip" 
can be resolved and integrated into the fluid-structure interaction dynamic problem as 
external loads. 

Thus, the motions of a large structure on or near a free surface (where viscosity 
effects can be ignored) could be obtained using the complete potential flow theory 



34 Le\vis. E.V. (ed.). Principles of Na\'al Architecture (l ot. Ill - Motions in Ifaves and Controllability). 
Second Revision. SNAME. Jersey City . NJ. 1989. 

35 Salveson. N.. Tuck. E.O.. and Faltinsen. 0.. "Ship Motions and Sea Loads". SKAXIE Transactions. Vol. 
78. 1970. pp. 250-279. 
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(linearized or fully nonlinear). It computes the forces exerted on the structure by the 
incident pressure (the external "loading" force - the bubble pulse for a whipping ship 
subjected to an underwater explosion bubble), the diffraction pressure (scattering of the 
incident pressure from the "rigid" structure), and the radiation pressure (wave radiation 
from the structure as it moves in each of its degrees of freedom) Thus, the total velocity 
potential (O) at any point in the fluid in the presence of a moving structure (either three- 
dimensional if panel methods are used, or two-dimensional if strip theory is used) can be 
written as a superposition of the velocity potentials due to incident, diffracted, and 
radiated pressures as follows: 



where <t>° is the velocity potential due to the incident pressure, O* is the velocity potential 
due to the diffracted or scattered pressure from the /th panel (where M is the total number 
of panels making up the total surface S), and <tT is the velocity potential due to the 

radiated pressure caused by motion of the /th panel in the yth direction (where N is the 
number of degrees of freedom for each panel) . 




M M N 



(2-54) 



1=1 1=1 j=i 
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Figure 2-17 Two-dimensional strip theory 
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For most marine applications, rigid body translations are ignored and the solution 
of the equations of motion is carried out in the frequency domain, assuming harmonic 
motion of the structure centered around an equilibrium position The incident and 
diffracted velocity potentials are computed on the assumption that the structure is fixed in 
its equilibrium position, subjected only to the given external loading (usually wave 
loading). The radiated velocity potentials are then computed assuming that the structure 
undergoes one simple harmonic motion at a time (i e at each frequency in each degree of 
freedom) in calm water (see Chakrabarti 36 , Faltinsen 37 , for example). For assumed simple 
harmonic motion of the form 

x = Xe“' (2-55) 

where x is the displacement vector, X is the displacement magnitude vector, and © is the 
frequency, the total hydrodynamic force on the structure (calculated with potential flow 
theory) has the form: 

F = -to : X M mn + icoX C mn (2-56) 

in which F n is the force on the structure in the //th degree of freedom, and X n is the 
displacement of the structure in the //th degree of freedom. is called the added mass 
matrix (hydrodynamic added mass in the //th degree of freedom due to unit displacement 
in the //7th degree of freedom), and thus the first term on the right hand side represents the 
hydrodynamic forces which are in phase with acceleration. C m is called the radiation 

damping matrix (hydrodynamic damping in the //th degree of freedom due to unit 
displacement in the ///th degree of freedom), and thus the second term on the right hand 
side represents the hydrodynamic forces which are in phase with velocity. 

For applications where simple harmonic motion about a fixed equilibrium position 
is not a reasonable assumption (e.g. when rigid body translations are important), the 
separate accounting for added mass and radiation damping forces as above is normally not 
possible (except in the case of steady rigid body translation). For cases with unsteady 
rigid body translation in addition to harmonic vibration (e.g. whipping of a submerged 
submarine subjected to an underwater explosion bubble), the calculation of the total 
hydrodynamic forces must be accomplished in the time-domain. Application of potential 
flow theory in the time-domain to predict general motions of structures is not new and 
fundamentally straight forward, but has had limited application for complex geometries 
and loadings because of limited computational resources. Discussions of direct time- 



36 Chakrabarti. op. cit 
37 Faltinsen. op. cit. 
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domain solution techniques are presented by various authors for various applications 
(Beck and Magee 38 , Ogilvie 39 , for example) 

2.4.2 Viscous fluid damping. 

As mentioned previously, dynamic energy can be lost into the surrounding fluid 
medium through forces associated with fluid viscosity. First, viscous shear forces occur 
within the viscous boundary layer surrounding the moving structure (so-called "skin 
friction" drag). Second, separation of the boundary layer from the structure leads to 
additional forces on the moving structure (so-called "form" drag). Together, boundary 
layer shear forces and boundary layer separation forces constitute the viscous fluid forces 

The extent to which each of these forces is important is a function of the geometry 
of the body, and the relative fluid motion around the body (characteristic amplitude of 
relative fluid motion and characteristic frequency). Where geometry and fluid flow around 
the structure are such that the boundary layer remains attached (e.g. streamlined bodies at 
low angles of attack), boundary layer shear forces tend to dominate. However, for "bluff' 
bodies such as circular cylinders and "flat" plates, boundary layer separation readily occurs 
which tends to produce forces greatly exceeding the boundary layer shear forces. 
Numerous experiments have been conducted on various geometries giving indication 
under which conditions each of the forces would be significant, and estimation of the 
forces. 

The largest body of data exists for harmonic transverse oscillation of horizontal 
cylinders for which two dimensional flow around the body can be assumed. As discussed 
previously, flow around a two-dimensional body is typically characterized by several non- 
dimensional parameters, including the Keulegan-Carpenter number and the Reynolds' 
number (or alternatively the viscous-frequency parameter, (3). Anaturk 40 and Anaturk et. 
al. 41 present regions in the KC - P plane where two dimensional flows possess different 
characteristics For example, when KC < 1 and P < 10’ — 1 0 4 (Rn < 1 O’ — 10 4 ), the 
boundary layer flow typically remains attached for the case of a smooth circular cylinder. 



38 Beck. R.F., and Magee, A.R.. "Time-domain Analysis for Predicting Ship Motions". Dynamics of 
Marine Vehicles and Structures in lf'a\’es. Elsevier Science Publishers. B.V.. 1991. pp. 49-64. 

39 0gilvie. T F.. "Recent Progress Toward the Understanding and Prediction of Ship Motions". 
Proceedings of the Fifth Symposium on Naval Hydrodynamics. Office of Naval Research. Washington. 
DC. 1964. pp 3-128. 

40 Anaturk. A.R.. "An experimental investigation to measure hydrodynamic forces at small amplitudes and 
high frequencies”. Applied Ocean Research. 13. 4, 1991. pp. 200-208. 

41 Anaturk. A.R.. Tromans. P.S., van Hazendonk. H.C.. Sluis. C M., and Otter. A.. "Drag forces on 
cylinders oscillating at small amplitude: a new model". Journal of Offshore Mechanics and Arctic 
Engineering. 114. 1992. pp. 91-103. 
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As KC and (3 (or Rn) increase, the turbulence in the boundary layer increases and flow 
separates. 

For the case where the two dimensional flow does not separate, the flow (and 
resulting forces) can be satisfactorily predicted using an attached laminar boundary layer 
theory (Newman 42 , Stuart 43 , Batchelor 44 , for example). This theory predicts the drag 
forces resulting from the combined normal and tangential stresses in the fluid using 
Navier-Stokes equations. 

For a body o f general shape, the forces associated with separated flow, and the 
point at which transition to separated flow occurs (i .e. the stability of the boundary layer) 
is dependent upon many factors. Many attempts have been made to solve separated flow 
around marine structures numerically (generally limited to two dimensional flow). 
Examples of various "state-of-the-art" methods include: 

(a) Single vortex method (Faltinsen and Sortland 45 ) 

(b) Vortex sheet model (Faltinsen and Pettersen 46 ) 

(c) Discrete vortex method (Sarpkaya 47 , Clements 48 ) 

(d) Vortex-in-cell method (Stansby and Dixon 49 ) 

(e) Combination ofChorin's method and vortex-in-cell method (Chorin 50 , 
Smith and Stansby 51 ) 

(f) Navier-Stokes solvers (Lecointe and Piquet 52 ) 

Additionally, attempts have been made to include vortex shedding models in numerical 
solutions for potential flows and thus combine wave radiation with the effects of viscosity 



42 Newman, op. at. 

43 Stuart. J.T., "Double boundary layers in oscillatory viscous flow ", Journal of Fluid Mechanics. 24. 

1966, pp. 673-6X7. 

44 Batchelor. G.K... . hi Introduction to Fluid Dynamics. Cambridge University Press. London. 1973. 
45 Faltinsen. O.M.. and Sortland, B.. "Slow-dnft eddy making damping of a ship". Applied Ocean 
Research. 9. 1. 19X7. pp. 37-46. 

46 Faltinsen. O.M., and Pettersen. B., "Application of a vortex tracking method to separated flow around 
marine structures". Journal of Fluids and Structures. 1, 19X7. pp. 217-237. 

47 Sarpkaya. T., "Computational Methods With Vortices - The 19XX Freeman Scholar Lecture". Journal of 
Fluids Engineering. Vol. 1 1 1. March 19X9. pp. 5-45. 

48 Clements. R.R.. and Maull, D.J.. "The representation of sheets of vorticity by discrete vortices". 

Progress in Aerospace Science. Vol. 16. No. 2. 1975. pp. 129-146. 

49 Stansby, P.K.. and Dixon. A.G., "Simulation of flow s around cy linders by a Lagrangian vortex scheme". 
Applied Ocean Research. Vol. 5, No. 3. 19X3. pp. 167-17X. 

50 Chonn. A.J.. "Numerical study of slightly viscous flow". Journal of Fluid Mechanics. 57. 1973. 
pp. 7X5-796. 

5l Smith. P.A.. and Stansby. P.K.. "Impulsively started flow around a circular cylinder by the vortex 
method". Journal of Fluid Mechanics. 194. 19XX. pp. 45-77. 

52 Lecointc. Y.. and Piquet. J. "Compact finite-difference methods for solving incompressible 
Navier-Stokes equations around oscillatory bodies". I 'on Harman Institute for Fluid Dynamics Lecture 
Series 1985-04, Computational Fluid Dynamics. 19X5. 
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(Yeung and Wu 53 , Wong and Calisal 54 , for example). While these methods have 
documented satisfactory agreement with experiment in some cases, they are generally not 
considered robust enough (numerical instabilities, etc.) to be applied with complete 
confidence to new problems where there is no guidance from experiments 55 . Thus, 
numerical methods to characterize separated flow are often used in conjunction with 
experimental methods. 

Morison's equations In practice, when drag forces are important (and wave 
radiation may be ignored), and two-dimensional oscillatory flow can be assumed, a 
Morison's equation or modified Morison's equation can be used to predict hydrodynamic 
forces . Originally proposed by Morison et. al 56 for the prediction of wave forces on 
vertical piles which extend from the bottom through the free surface, the Morison 
equation has been modified to include application to cylinders of varying orientation under 
various loading conditions, including cases where the cylinders are free to move under the 
imposed fluid loading. 

Morison's equation for a stationary ’ cylinder in an oscillating flow field includes 
the effect of fluid inertia forces and assumes drag forces to be a quadratic function of 
velocity: 

F = pAC m u + ^pDC d u|u| (2-57) 

where F is the hydrodynamic force (per unit length of cylinder - similar to "strip theory", 
the total force on the cylinder can be found by integrating over the length of the cylinder), 
A is cross sectional area, D is diameter, u is fluid velocity ( u is fluid acceleration), 
coefficient C m is called the inertia coefficient (the first term being proportional to 
acceleration is equivalent to an inertia force), and coefficient C d is called the drag 
coefficient (the second term being proportional to velocity (squared) is equivalent to a 
drag force). Figure 2-18 illustrates the parameters used in Morison's equation. 

Morison's equation for an oscillating cylinder in a stationary’ fluid includes the 
effect of fluid added mass forces and again assumes drag forces to be a quadratic function 
of velocity, giving the reaction forces on the cylinder: 

53 Yeung. R.W.. and Wu, C.. "Viscosity effects on the radiation hydrodynamics of two-dimensional 
cylinders". Proceedings of the l Oth International Conference on Offshore Mechanics and Arctic 
Engineering. ASME. 1991. pp. 309-316. 

54 Wong, L.H.. and Calisal. S.M.. "A numerical solution for potential flows including the effects of vortex 
shedding". Proceedings of the 1 1th International Conference on Offshore Mechanics and Arctic 
Engineering. ASME. 1992, pp. 363-368. 

55 Faltinsen. op. cit. 

56 Morison. JR. O'Brien. M.P.. Johnson. J.W.. and Schaaf. S.A.. "The force exerted by surface waves on 
piles". Petroleum Transactions. A1ME. 189, 1950. pp. 149-157. 
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(2-58) 



F 1 = -pAC A x-^pDC d x|x 

where x and x are the acceleration and velocity of the cylinder, C A is called the added 
mass coefficient, and Cj is the drag coefficient for and oscillating cylinder in still water. 




Figure 2- 1 8. Parameters used in Morison's equation 

Coefficients C m , C d , C A , and C d are normally experimentally determined for 
various geometries and loading conditions (usually as functions of KC and P or Rn), and 
have been determined in a large number of experimental studies (Bearman et. al. 57 , 
Sarpkaya 58 , Chakrabarti and Cotter 59 , for example). 

When a cylinder is free to move in a moving fluid, Morison's equation can be 
written in an extended form by combining equations (2-57) and (2-58): 

F = pAC m ii + ^pDC d u|u|-pAC A x-^pDC d x|x| (2-59) 

This form is known as the independent flow fields model 60 and is obtained as the 
superposition of the two independent flow fields; a far field due to the water particle 
motion (i.e. the external loading) and relatively unaffected by the presence and motion of 

57 Bearman. P.W.. Lin. X.W.. and Mackvvood. P.R.. "Measurement and prediction of response of circular 
cylinders in oscillating flow”. Proceedings of the 6th International Conference on the Behavior of 
Offshore Structures (BOSS), 1992. pp. 297-307. 

58 Sarpkaya. T.. "Force on a cylinder in viscous oscillator) flow at low Keulegan-Carpenter numbers". 
Journal of Fluid Mechanics, Vol. 165. 1986. pp. 61-71. 

59 Chakrabarti. S.K.. and Cotter. DC,, "Hydrodynamic coefficients of mooring tower". Proceedings of the 
15th Annual Offshore Technology } Conference , OTC 4496 , 1983. pp. 449-458. 

60 Chakrabarti, op. cit .. p. 187. 
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the structure, and a near field resulting from the motion of the structure Coefficients C w 
and Cj may be obtained from experiments where the cylinder is held fixed in oscillating 
fluid, and coefficients C A , and Cj may be obtained from experiments of an oscillating 
cylinder in otherwise still fluid Alternatively, the coefficients could be determined from 
empirical relations (using KC and Rn, for example): 



C C = f 



KC, =^,Rn, 

D 1 




(2-60) 



c c = f 

'-'A d 1 



KC„ = 



x,T 



x ,D 

, Rn =^— 



(2-61) 



D v 

where u 0 ,x 0 are the amplitudes of water particle and structure velocity, respectively, 

T,T‘ are periods of oscillation of water particle and structure, respectively, and subscripts 
/ and n refer to the far field and near field, respectively. 

An alternative to the independent flow fields model is the relative velocity model , 
in which forces are written in terms of the relative motion between the structure and the 
fluid, in which single coefficients are assumed to apply: 



F = P A C m (u- >0 + pAx + :|-pDC d (u-x)|u-x| (2-62) 

In the relative velocity model, the coefficients may be determined using empirical relations 
(using KC and Rn, for example): 



c m ,c d = f 






KC = Vr ° Tr 



,Rn.=^ 



J 



(2-63) 



D v 

where subscript r refers to the relative motion, v r is the relative velocity between the 
structure and the fluid (v r0 is the amplitude of the relative velocity), and T r is the 
combined period of the relative motion. The relative velocity model has seen extensive 
use in the analysis of hydrodynamic forces on drag dominated offshore structures (Laya et 
al. 61 , Spidsoe and Karunakaran 62 , for example). 

Thus, the Morison equation can be used under various conditions to predict the 
hydrodynamic forces on structures (or portions of structures) where drag forces and fluid 



61 Laya. E.J.. Conner, M.. and Sunder. S.S.. "Hydrodynamic forces on flexible offshore structures". 
Journal of Engineering Mechanics, ASCE. March 1984. pp. 433-448. 

62 Spidsoe. N.. and Karunkaran. D.. "Nonlinear effects of damping to dynamic amplification factors for 
drag-dominated offshore platforms". Proceedings of the 1 1th International Conference on Offshore 
Mechanics and Arctic Engineering, ASME. 1992. pp. 231-239. 
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inerm forces are important An appl, cation of Morison's equation to the development of a 
model tor submarine whipping will be made in chapter 3 



47 



3. MODEL FOR SUBMARINE WHIPPING INCLUDING DAMPING EFFECTS 



3.1 Introduction. 

In order to gain some quantitative insight into the relative amount of damping 
which could be provided by each of the mechanisms discussed in the previous chapter, it is 
first necessary to develop & general model for submarine whipping, which may include any 
or all of the specific mechanisms One way of doing this would be to develop a simplified 
model of each of the mechanisms based upon the fundamental principles of mechanics , 
some of which were used in the discussion in the previous chapter, and utilize a flexible 
formulation for the equations of motion which permits inclusion of selected mechanisms. 
Using this approach, any individual mechanism could be evaluated separately for its 
expected effect on total damping 

The equations of motion. The dynamic behavior of any structure is described by 
the equations of motion, which are simply extensions of Newton's second law of motion. 
Fundamentally, an equation of motion can be written as a sum of forces which resist 
motion of the structure (i e. the inertial forces (F incmal ), the damping or dissipative forces 

(FdampmgX and the elastic or stiffness forces (F elasuc )), and forces which induce motion (i.e. 
the external or loading forces (F loading )). Thus, an equation of motion may be written in the 
general form: 

F + F +F=F fl-H 

inertial damping elastic loading \ / 

Inertial forces are those forces through which kinetic energy is stored, and are 
dependent upon acceleration and a mass coefficient (F inemal (M, x)). Elastic forces are 
those forces through which potential energy is stored and are generally dependent upon 
displacement and some elastic or stiffness coefficient (F elastlc (K,x)). As discussed 
previously, damping forces are those forces related to the dissipation of energy, and 
depend upon the mechanisms through which energy is dissipated, generally dependent 
upon velocity and some damping coefficient (F damping (C,x)). It should be noted that, in 

general, these dependencies are neither constant nor linear 

The external loading forces include all forces which are "externally" applied to the 
structure. These forces may themselves be thought of as being inertial, dissipative, or 
elastic in nature. Thus, equation (3-1) may be rewritten by redefining the inertial, 
damping, and elastic terms such that they also include the "external" loading forces. For 
example, for a system consisting of a structure moving in a fluid, equation (3-1 ) may be 
rewritten 
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F +F 1 +F 1 — o 

inertial damping 1 elastic ^ 

where F.len.ai represents all inertial forces (inertia of the structure and inertia of the 



(3-2) 



surrounding fluid), represents all damping or dissipative forces (energy dissipation 

within the structure and energy lost into the surrounding fluid), and Fj lasti represents all 
elastic forces (structural stiffness and hydrodynamic buoyancy forces) 

In general, the terms associated with inertial forces, damping forces, and elastic 
forces may be written in any order as long as all forces are accounted for (and as long as a 
proper sign convention is maintained) In keeping with the classification of damping 
forces discussed in chapter 2, the equations of motion for a whipping ship might therefore 
be written in the form 



where; 



^hull ■*" F 



int emal 



+ External ~ 0 



(3-3) 



F — i F + F 4 - f \ 

r hull l r inertial ^damping ^elastic J 

F = ( f 4- F 4- F ) 

internal l inertial A damping elastic J 

F = (f + F +F ) 

external l inertial damping elastic J 



hull 



int emal 



external 



(3-4) 



where F exlemal would include all external forces on the hull Thus, the dynamics of a 
whipping ship may be thought of as a (coupled) superposition of the dynamics of the hull, 
any internal components, and the external fluid medium. Equations (3-3) and (3-4) 
represent one form of the complete, coupled, nonlinear equations of motion for a 
whipping ship. 

Using this approach to the equations of motion, each term may be represented by 
an independent formulation (which may involve specific assumptions of linearity for that 
term), so long as the terms remain coupled via a single solution technique. As an example, 
the dynamic behavior of a slender ship hull (such as a modem submarine) might be 
calculated based upon a linear simple beam finite element model, while the external 
(hydrodynamic) forces might be calculated based upon a nonlinear three-dimensional panel 
method. As will be discussed subsequently, even if the hull itself can be assumed to 
respond linearly, the nonlinearities associated with the external (hydrodynamic) forces 
require the solution of the equations of motion in the time-domain It is the existence of 
the coupling of the terms within the equations of motion which require the solution of the 
combined equation by a single numerical algorithm. 
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3.2 Formulation of the dynamic analysis procedure. 

In general structural dynamics, the selection of the most appropriate method to 
solve the equations of motion depends upon the type of loading on the structure, the 
extent to which the equations of motion can be uncoupled and linearized, and the 
existence of rigid body motions Nonlinearities which must be addressed include 
geometric nonlinearities (large rotations), material nonlinearities (nonlinear material stress- 
strain relations), nonlinearities associated with the external loading (hydrodynamic 
nonlinearities), nonlinear boundary conditions, and so on. The extent to which 
nonlinearities can be "linearized" depends upon the assumptions made in developing the 
physical model for the motions of the particular structure. For ship whipping, even if 
material and geometric nonlinearities can be effectively linearized (by assuming linear 
elastic stress-strain relations and small rotations - reasonable assumptions for the case of 
elastic whipping of a steel ship hull), hydrodynamic forces can be expected to be highly 
non-linear (as discussed in chapter 2), and thus complete linearization of the equations of 
motion is generally not possible Additionally, for the general dynamic response of a 
whipping ship subjected to an underwater explosion bubble, it could be expected that there 
would be a significant potential for rigid body motion. For these reasons, the equations of 
motion must be solved directly in the time-domain (utilizing a time-stepping numerical 
algorithm). 

Implementation of time-domain analysis procedures for the dynamic response of a 
structure in a fluid medium is discussed in great detail in a large number of texts and 
journals articles. Bathe 63 , Clough and Penzien 64 , and Argyris and Mlejnek 65 provide 
particularly thorough discussions of the solution of nonlinear equations of motion by direct 
time integration methods. Direct time integration means that the equations of motion 
would be solved using a numerical step-by-step procedure. The term "direct" generally 
implies that there is no transformation of the equations into a different form prior to 
numerical integration, although using the formulation of the equations of motion discussed 
above, it is feasible to utilize a coordinate transformation, so long as it is conducted at 
each time step (this will be discussed in greater detail later). In essence, instead of trying 
to satisfy equilibrium at only one distinct time / (as for static equilibrium), direct time 
integration aims to satisfy dynamic equilibrium at all discrete time intervals At apart This 
implies that, basically, dynamic equilibrium (which includes the effects of inertia, damping 



63 Bathe. op. cit. 

64 Clough and Penzien. op. cit. 

65 Argyris. J.. and Mlejnek. J.P.. Dynamics of Structures. Elsevier Science Publishers. The Netherlands. 
1991 . 
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and elastic forces) is sought at discrete time points within the interval of solution 
Typically, for direct time integration algorithms, the solution at the "next" required time, 
t+At, is calculated based upon the solutions at all previous times, and assuming an 
incremental variation in displacements, velocities, and accelerations. It is the form of the 
assumption on the variation of displacements, velocities, and accelerations within each 
time interval that determines the accuracy, stability, and required time of the solution 
procedure 

For direct time integration methods, the incremental equation of motion can be 
written as an extension of equation (3-2) by subtracting the equilibrium equation at time t 
from that at time t+At: 

fpu-At pt j,(pt+At _ pi j i ( p<^‘ _ c' )_n 

L inertial inertial J l * damping A damping J l elastic elastic J v*' 

As discussed above, inertial forces are dependent upon acceleration and a mass coefficient 
( F.nemai (Mi x)), damping forces are dependent upon velocity and some damping coefficient 

(^damping (C,x)), and elastic forces are dependent upon displacement and some elastic or 
stiffness coefficient (F elas[ic (K,x)). Thus, if coefficients M, C, and K can be related to 
acceleration, velocity and displacement (respectively), then equation (3-5) contains only 
three unknowns at time t+At (the solution at time t would be known) and may be solved if 
two of the three are known as relations to the third. 

The Newmark time integration method. The Newmark integration method 66 , 
which is one of the most reliable direct time integration schemes, and whose stability 
condition is well known, can be employed quite appropriately for this application (see also 
Bathe 67 , Argyris and Mlejnek 68 , and Doyle 69 ). In the Newmark method (often referred to 
as the Newmark (3 Method), iterative calculation is carried out at each time step and is 
moved to the next time step after an equilibrium condition is reached. 

The basic equations of the Newmark method are derived by assuming Taylor 
expansions for displacements (x) and velocities (x ): 

x, + a. =x,+x l At + x t ^- + p(x l+Al -x t )At 2 (3-6) 

x l+ At = x, + x,At + y(x t+At - x t )At (3-7) 



66 Nevvmark. N.M., M A method of computation for structural dynamics." American Society' of Civil 
Engineers Transactions. Vol. 127. 1962, pp. 601-630. 

67 Bathe. op. cit. 

68 Argv ris and Mlejnek. op. cit.. pp. 485-488. 

69 Doyle, J.F., Static and Dynamic Analysis of Structures. Kluwer Academic Publishers. Boston. 1991. pp. 
353-357. 
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where P and y are coefficients to evaluate contributions of remainders. Thus, the 
displacements and velocities are approximated based upon "assumed" accelerations at 
each time interval (thus requiring an iterative approach), and based upon the solution at 
the previous time interval. Newmark showed that y and (3 are not entirely arbitrary and 
their selection determines the stability and accuracy of the numerical solution Newmark 
found that unless y is greater than or equal to 0.5, spurious damping is introduced into the 
equations of motion. Additionally, unless (3 is greater than or equal to 0.25, there is an 
incorrect maximum velocity response. Further, with y = 0.5 and P = 0.25, the scheme 
leads to unconditionally stable solutions 

To implement the Newmark method in the analysis of the dynamics of a structure 
subject to nonlinear loading forces (such as a whipping submarine), the equations of 
motion are first put in the form: 

^resisting ^loading (-^“^) 

where F resisting refers to all those forces associated with the structure itself which resist 
motion (i.e. F hull as defined in section 3.1) and F, oading refers to all (other) forces which 
induce motion (i.e. F mtemal + F extemal as defined in section 3.1). Defining F resisting for a 
structure as: 

=Mx + Cx + Kx (3-9) 

the equations of motion for a structure can be written in the general form 

Mx + Cx + Kx = F loadmg (3-10) 

or in incremental form: 

M(x tiit - x t ) + C(x t ^ t - x, ) + K(x t ^ t - x t ) = F loaJmgt+At - F loaJmg[ (3-11) 

From equation (3-10), the acceleration required for dynamic equilibrium to exist at 
time t+At can be written: 

-{Cx,. 4 , +Kx,.j} (3-12) 

Because F loading is generally also dependent upon the accelerations (and displacements and 

velocities by equations (3-6) and (3-7)), it is evident that an iterative procedure must be 
implemented at each time t+At in order to obtain equilibrium by equation (3-12). 
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Newmark 70 proposed a general method for the solution of structural dynamics 
which may be summarized as follows: 

1. Assume initial values of displacement, velocity, and acceleration (i.e. at 

each node) for time t = 0. 

2 For each time step: 

(a) Assume values of acceleration for the current time step. 

(b) Compute displacement and velocity for the current acceleration 
(from equations (3-6) and (3-7). 

(c) For the current displacement, velocity, and acceleration, compute 
the loading and resisting forces in accordance with the physical 
model for the particular application 

(d) Compute the acceleration required to obtain dynamic equilibrium 
with these loading and resisting forces (from equation (3-1 1)). 

(e) Compare the required acceleration with the assumed acceleration 
for the current time step If these are the same (i.e. within 
tolerance), the calculation for this time step is complete. If these 
are different, repeat the calculation (i.e. perform another iteration 
through steps (a) - (e)) with a different value of assumed 
acceleration (the derived acceleration could be used as the new 
assumed acceleration). 

3 Move onto the next time step (t+At), assuming the final acceleration for the 

last time step (t). 

It should be briefly noted here that the modeling of internal damping devices such 
as structures or liquids requires only the addition of N equations of motion to the set 
defined above (for an internal damping device with N participating modes of response). 
Thus, in addition to the requirement for equations (3-12) to be satisfied at each time t+At, 
the following equations must also be satisfied at each time t+At: 



where subscript n is for each mode of the internal damping device (up to mode N). Thus, 
if there are M degrees of freedom required to represent the displacements of the hull, and 
N (modal) degrees of freedom required to represent the (modal) displacements of the 
internal device, then a total of M^N simultaneous equations must be satisfied using the 
form of equations (3-12) and (3-13). 

70 Ne\vmark. op at.. pp. 606-61 1. 
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To implement the Nevvmark method in the analysis of the dynamics of a structure 
subject to loading that is highly nonlinear (such as a whipping submarine subject to 
nonlinear hydrodynamic forces), a slight modification to the above method must be made 
Because of the highly nonlinear loading forces, the equations of motion must be solved at 
each time step using an iterative process (such as a Newton-Raphson process) which 
"guarantees" convergence. Without such an iterative process, the solution may not 
converge and equilibrium may not be (numerically) obtained 

Suzuki and Yoshida 71 discuss application of the Nevvmark method with solution at 
each time step by an iterative Newton-Raphson method In the case of elastic ship 
whipping with nonlinear hydrodynamic loading, the specific formulation of the Newton 
iteration may proceed by using equations (3-6), (3-7) and (3-1 1), and expressing 
acceleration and velocity at time t+At in terms of displacement at time t+At 



L l-At 



— x, + 



i 



(3-14) 



"l-At 



JL( X ^- Xi)+ ,_X 

pAt lAt ' l p 






x, + 



1-— x At 

23 J 



(3-15) 



In order for equilibrium to exist at time t+At, equations (3-1 1), (3-14) and (3-15) must be 
simultaneously satisfied Residuals (which must vanish for equilibrium to exist) may be 
defined for equations (3-1 1), (3-14) and (3-15): 
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The Newton iteration is carried out such that the unsatisfied equilibrium at iteration j is to 
be satisfied at iteration j+1. In other words, all residuals at iteration j+1 must 
simultaneously vanish To obtain expressions for residuals at iteration step j+1, the 



71 Suzuki, H.. and Yoshida. K.. "Three dimensional nonlinear dynamic analysis method of underwater line 
stmcture and its validation". Proceedings of the l Oth International Conference on Offshore Mechanics 
and Arctic Engineering (OMAF.). ASME. 1991. pp 87-94. 
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residuals at iteration step j may be expanded in linearized Taylor expansion (in matrix 
form): 
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displacement, velocity, and acceleration corrections to be made for the j+1 th iteration at 
time step t+At. An expression for the displacement correction (5x) may be obtained by 
substituting the last two equations in (3-19) into the first, giving: 
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(3-20) 

Solving equation (3-20) for 5x and substituting this term back into equation (3-19) gives 
expressions for improved displacement, velocity and acceleration (i.e. displacement, 
velocity and acceleration for iteration step j~ / based upon iteration step j) 
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The Jacobian matrices, 
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numerically differentiating the loading vector. 



, may be approximated at each iteration step by 
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With these modifications, the calculation of dynamic equilibrium using the 
Newmark method for the case of highly nonlinear loading could proceed as follows 
1 Assume initial values of displacement, velocity and acceleration (i.e 
vectors of nodal values). 

2. For each time step: 

(a) Assume initial values of displacement, velocity, and acceleration 
equal to those of the last time step 

(b) Calculate the loading forces for the current displacement, velocity, 
and acceleration using the chosen physical model. 

(c) Calculate the Jacobian matrices (numerically) for the current 
loading forces. 

(d) Calculate the residuals using equations (3-16) through (3-18). 

(e) Solve for the displacement correction (5x) using equation (3-20) 

(0 Calculate improved values of displacement, velocity, and 

acceleration using equations (3-21) 

(g) Reiterate through steps (b) through (g) until the residuals are within 
a specified tolerance (i.e. a vector norm tolerance). 

Move onto the next time step, assuming the final displacements, velocities, 
and accelerations of the previous time step 
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3.3 



Definition of dynamic forces. 



3.3.1 Hull structure forces. 

As discussed in section 3.1, the forces associated with the dynamics of the hull 
structure may be represented by an independent formulation from that used to represent 
the internal and external forces. This formulation may involve specific assumptions of 
linearity which may be applied only to the dynamics of the hull structure 

3.3. 1.1 Finite element beam model. 

In representing the whipping response of most conventional ships to low frequency 
bubble pulse loading, it has been found that a simple beam model of the hull structure has 
provided good results 72 . This idealization is supported by recognizing that most 
conventional ships (particularly submarines) have lengtlVbeam ratios greater than eight. 
Thus, the low frequency vibrational modes occur as beam-like flexural motions, which are 
represented well by simple beam theory. It should be noted, however, that nonlinear hull 
response (i e. where the assumptions of linear elastic stress-strain relations and small 
rotations cannot be made), makes the utility of the simple beam theory limited For 
nonlinear hull response, a fully nonlinear, three-dimensional finite element technique would 
be more appropriate. Therefore, for this analysis procedure, the assumptions of linear 
elastic stress-strain relations and small rotations will be made in order to justify the use of 
the simple beam model. 

Because sectional properties vary gradually along the length of the hull, a finite 
element beam model could be used to represent the structural properties of the submarine 
(figure 3-1) The model consists of a specified number of beam elements (equal to the 
number of hull stations specified on the ship's drawings). Most ships have 19 stations 
specified, and thus would be modeled with 19 beam elements. A three-dimensional beam 
model will be used because of the directionality of the loading of the underwater 
explosion, as well as the lack of complete axial symmetry of submarines (i .e properties 
differ in each plane of symmetry). 

Most methods for structural finite element analysis develop matrices to represent 
the stiffness and mass properties of each finite element , and then combine element matrices 
into global or structure matrices. Modeling of beams using a finite element technique, and 
assignment of the properties of each finite element is discussed in numerous sources 



72 Hicks. op. cit., p. 395. 
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(Smith 75 , Bathe 74 , Doyle 75 , Hughes 76 , for example). An appropriate finite element analysis 
procedure which could be used to approximate the hull forces associated with submarine 
hull whipping could be outlined as follows: 

1 Define the model (coordinate system, nodes, element types, etc.). 

2 Determine each element's stiffness and mass matrices ( k \ m f ) and 
transform them to hull (or global) coordinates. 

3 . Assemble the hull or global stiffness and mass matrices ( K , M ) by 
superposition of the element matrices. 

Approximate a reasonable hull damping matrix (C) by applying appropriate 
Rayleigh coefficients, as discussed in sections 2.2.3 and 3.3. 1.2 
(C = aM+pK). 

Combine the hull mass, damping, and stiffness matrices to obtain the total 
hull force as described in section 3.1 (F hull = Mx + Cx + Kx). 

This procedure will be incorporated into an overall computational algorithm and used for 
analyses in the next chapter. 



4. 



5. 





Elastic bean Lunpeoi 

elements masses/nodes 



Figure 3-1. Structural model for submarine hull whipping. 

Figure 3-2 illustrates a two-dimensional, two-node beam element with vertical 
translation and rotation as the nodal displacements. It should be noted that because a 
submarine should be classified as a relatively thick beam (i.e. the cross-sectional 
dimensions are not small compared to the length), the effects of rotary inertia and shear 
deformation at each element should be accounted for in the finite element model. 



7, Snwth. op. cit. 

?4 Bathe. op. cit. 

75 Doyle, op. cit. 

76 Hughes, O F.. Ship Structural Design , SNAME. Jersey City. NJ. 1988. 



58 



To be rigorous, a three-dimensional beam element could be used which also considers 
additional nodal displacements in the horizontal plane as well as torsion and axial 
displacements (a total of 6 degrees of freedom at each node). Figure 3-3 illustrates a 
general three-dimensional, two-node beam finite element with 6 degrees of freedom at 
each node. It should be noted that axial displacements can generally be ignored for 
submarine whipping due to the lack of axial loading provided by the surrounding fluid. 
Additionally, torsional displacements may also be ignored, although it could be argued that 
torsional loading could become important under some conditions in accounting for 
potential asymmetry associated with the underwater explosion bubble loading on external 
hull appendages such as the sail and planes. Therefore, axial and torsional degrees of 
freedom will be omitted for this analysis, although the entire 3-dimensional beam element 
will be used for illustration in the following discussion. 



y 




Figure 3-2. Two-dimensional, two-node beam finite element. 




Figure 3-3. General three-dimensional, two-node beam finite element. 
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In order to simplify the assembly of the global equations of motion for the hull, the 
equations of motion for each beam finite element may be written in the partitioned form 



M 
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M 21 
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where subscripts I and 2 refer to nodes 1 and 2 of the element, and each partition is a 6x6 
submatrix or 6 element subvector For the general three-dimensional beam element, the 
subvectors for nodal displacement and nodal force would be written in the form 
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where subscript / refers to the node ( 1 or 2), x is nodal translation, 0 is nodal rotation, f is 
the (external) force acting on the node, and m is the (external) moment acting on the node. 
It will be shown subsequently that formulation of submatrices for mass and stiffness can be 
readily accomplished However, the definition of element damping matrices as in equation 
(3-22) is generally not possible, thus a formulation of global hull damping will be given in 
the next section which will approximate the overall energy dissipation within the hull 
structure. 

The stiffness matrix. The derivation of stiffness properties for a finite element is 
a lengthy process and will not be discussed in detail here Barltrop and Adams 77 and 
Hughes 78 provide formulation for the element stiffness matrix of a general three- 
dimensional, two-noded beam finite element. The element stiffness matrix can be written 
in a node-oriented partitioned form consistent with equations (3-22) and (3-23). Each 
partitioned submatrix (6x6) has terms ij which represent the element elastic force in the 
direction of freedom / when a unit displacement is applied to freedom j with all other 
freedoms held at zero displacement The stiffness submatrices consistent with equation 
(3-22) can be defined as shown in figure 3-4. 

The mass matrix. The simplest method of assigning mass properties of the 
ship/beam model is to assume that the mass of each section is lumped or concentrated at 



77 Barltrop. N.D.P.. and Adams. A.J.. Dynamics of Fixed Marine Structures, Butterworth-Hcinemann. 
UK. 1991. Appendix D 
78 Hughes. op. cit.. pp. 212-217. 
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the nodal points (figure 3 - 1 ). This is called the lumped mass formulation In this way, the 
lumped masses fill the diagonal terms of the element mass matrix (corresponding to the 
individual nodal degrees of freedom), while the off-diagonal terms are all zero. The 
matrix (or submatrix) elements would therefore equal the mass of the node for elements 
corresponding to translational degrees of freedom, or the mass moment of inertia of the 
node for those elements corresponding to rotational degiees of freedom. The translation 
of a lumped mass accounts for translational inertia effects, while the rotation of the 
lumped mass accounts for rotary inertia effects. 

The lumped mass formulation has some very appealing properties. Firstly, it is 
easy to associate a physical model with the matrix (each degree of freedom at each node 
having its own mass). Secondly, the entire structure mass matrix is diagonal, leading to 
significantly fewer computations and storage requirements. The finite element mass 
submatrices as defined in equation (3-22), for a lumped mass beam finite element, could be 
written as shown in figure 3-5. 
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A is the cross-sectional area of the element, E is Young's Modulus, L is the length of the 
element, I z ,1 are the 2nd moments of area of the element cross-section about the z and 
y axes, G is the shear modulus. J is the polar moment of inertia of the element cross- 
section (about the x axis), and A w , A sz are the portions of the cross-sectional area over 
which the y and z-directed shear forces are assumed to act. 



Figure 3-4 Partitioned submatrices for stiffness of a general 3-D beam finite element 
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M 1; = M :1 = [o] (no off-diagonal terms for the lumped mass formulation) 
m, is the mass associated with node i 
I* >I y , .I a are the mass moments of inertia of node i about axes x,y,z. 

Figure 3-5. Partitioned submatrices for mass of a general 3-D beam finite element. 

It is also possible to take advantage of the finite element formulation to develop 
what is called a consistent mass formulation for a beam element. In contrast to the 
lumped mass matrix, the components of the consistent mass matrix cannot be associated 
with a physical model. More importantly, after assemblage, the global mass matrix will 
not be diagonal, but will include both diagonal and off-diagonal terms. For this reason, 
the consistent mass matrix formulation requires considerably more computational effort 
than the lumped mass formulation. While the consistent mass matrix formulation is 
generally superior in terms of minimizing errors in calculation of natural frequencies 79 , the 
errors associated with the lumped mass formulation significantly decrease as the number of 
nodes used to model the beam increases. Additionally, when external loads are being 
applied directly to the nodes (as will be discussed subsequently), the equations of motion 
behave better when the inertial forces are likewise lumped at the nodes. For these reasons, 
the lumped mass matrix formulation will be used in this submarine whipping model. 

3.3.1. 2 Hull damping model. 

Rayleigh damping. As mentioned previously, an approximation for the global 
damping matrix (C) associated with the hull structure can be made by applying 
appropriate Rayleigh coefficients to the global mass and stiffness matrices derived by 
superposition of the element mass and stiffness matrices: 



79 Doyle. op. cit . . p. 272. 
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C = otM + f3K 



(3-24) 



where M and K are the global mass and stiffness matrices, and coefficients a and (3 are 
the two damping proportionality constants. As discussed in section 2.2.3, coefficients a 
and (3 can be calculated if two modal damping ratios and £ b ) associated with two 
specific modal frequencies (a a and a b ) are known 



The downside to the Rayleigh (or proportional) damping matrix formulation is that 
the contribution of the higher frequency modes of hull flexural response (modal 
frequencies much greater than o b ) are effectively eliminated by the high damping imposed 



whipping, where only the lower few modes of flexural response are generally required to 
adequately describe the motion of the structure, this limitation can be considered 
acceptable as long as the higher frequency used for evaluation of the Rayleigh coefficients 
(o b ) is set among the higher modes expected to contribute significantly to the whipping 
response. Thus, the problem becomes one of calculating the modal (natural) frequencies 
for the lowest modes of hull whipping, applying experimentally determined (or reasonably 
approximated) modal damping ratios for two of the modes, and calculating the appropriate 
Rayleigh coefficients using equation (3-25). As discussed in section 2.2.3, the modal 
(natural) frequencies of a submarine structure can be calculated by solving the generalized 
eigenproblem for the undamped free vibration of the dynamic system: 



where o n is the natural frequency of vibration of the //th mode and <j) n is the mode shape 
vector of the mh mode. Data for experimentally determined (hull) modal damping ratios 
for the fundamental modes of submarine vibrations is lacking in the literature. However, 
some work has been done in the marine industry to estimate hysteretic and structural 
modal damping ratios for steel pile platforms for fundamental-mode flexural vibrations in 
random seas 80 , implying that hysteretic damping in thin steel cylindrical structures is 
typically less than 1% of critical (i.e. C, < 0.01) for the fundamental flexural mode. There 
is also some older data available in the literature for modal damping ratios of steel circular 
cylinder ship masts based upon modal "bump" testing 81 , suggesting that general steel 



80 Cook. M.F.. and Vandiver. J.K.. "Measured and Predicted Dynamic Response of a Single Pile Platform 
to Random Wave Excitation". Proceedings of the 1 4th Offshore Technology Conference. OTC 4285. 
1982. pp. 637-645. 

81 Smith, op. cit. 




(3-25) 



by the selection of the particular coefficients (see figure 2-4). In the case of elastic ship 




(3-26) 
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cylinders have (material and structural) modal damping ratios on the order of 2-3 % of 
critical for the lower flexural modes, although this data seems less reliable because of 
other unknown damping effects on the mast not considered (air and structural connection 
to the ship hull) Although submarines are not precisely the same as thin steel pile 
structures or ship masts, the general hysteretic damping mechanisms are basically the same 
and thus it would seem reasonable to extend this to data for a similar cylindrical steel 
structures to a steel submarine Thus, it would be reasonable to assume (as a first 
approximation) modal damping ratios (for the hu/f) on the order of 1% for the 
fundamental flexural modes of a steel submarine 

Thus, once the modal (natural) frequencies of the hull structure have been 
calculated by equations (3-26), reasonable Rayleigh coefficients could be calculated by 
equations (3-25) assuming modal damping ratios on the order of 0.01 Finally, the hull 
damping matrix to be applied in the finite element calculations can be calculated using the 
global mass and stiffness matrices and equation (3-24). 

This process will be incorporated into the computational algorithm and used in 
analyses in the next chapter. 

An alternative approach. "Dry mode" modal superposition. An alternative to 
developing an explicit damping matrix (C) using Rayleigh damping and solving the 
complete equations of motion, is to use the approximated hull modal damping ratios 
directly and a modal superposition approach As discussed in section 3.1, each vector 
term of the equations of motion may be represented by an independent formulation (which 
may involve specific assumptions of linearity for that term), so long as all the terms remain 
coupled via a single numerical solution algorithm. Although forces associated with 
external (fluid) loading may be nonlinear, and the equations of motion must therefore be 
solved in the time domain, it is still possible to solve for the forces associated with the hull 
structure using a modal superposition technique. This type of approach has been referred 
to as "dry mode" modal superposition 82 . In a "dry mode" modal superposition approach, 
all forces not associated with the hull are simply termed "loading" forces and thus carried 
on the right hand side of the equations of motion. 

{Mx + Cx + Kx}„ = F_*, (3-27) 

Thus, a coordinate transformation (of the generalized coordinates of the ship hull) may be 
conducted, so long as the transformation of the loading vector ( i.e. F mtemal and F extemal as 
discussed previously) is conducted at each time interval. 



82 Bishop. R.E.D.. and Price. W.G.. Hydroelasticity of Ships. Cambridge University Press. Cambridge. 
UK. 1979. 
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This formulation is consistent with the formulation of the Newmark time 
integration method presented in section 3.2 (i.e. equation (3-10)). Because the forces 
associated with the hull are to be modeled using assumptions of linearity (as discussed in 
the previous section), these forces may be redefined using the principles of modal 
superpostion outlined in section 2 2 3 The "dry" mode shapes are first found by solving 
the generalized eigenproblem for the undamped free vibration of the hull structure using 
equation (3-26). The modal equations of motion (equations (2-12)) are then written from 
equations (3-27) by applying the modal transformation, equations (2-10) and (2-13). The 
transformation of the modal force (i.e Q n in equations (2-13)) must be performed at each 
iteration at each time step of the direct time integration procedure, since the loading 
forces (F Ioad ) are generally nonlinear and coupled to the generalized displacements, 
velocities, and accelerations (i.e are functions of x,x,x). 

While this approach has some appealing qualities in terms of representation of the 
response of a whipping ship in terms of its fundamental flexural mode shapes, it affords 
little (if any) computational advantage over rigorous computation using the generalized 
coordinates in that two transformations must still be made at each iteration step of each 
time step. Thus, the fundamental benefit of modal decomposition is lost due to the 
requirement of maintaining the complex (i .e. nonlinear) external loading forces in terms of 
the generalized coordinate system. For this reason, this approach will not be pursued in 
this analysis and all computation will be carried out in terms of generalized coordinates. 

3.3.2 Internal forces. 

In order to analyze the potential contribution of discrete internal structures and 
sloshing liquids on the whipping response of a ship, it is necessary to assume some 
"reasonable" properties of the internal structure and liquid storage tanks as might be 
found on a "real" ship. Appropriate models to account for the dynamics of the internal 
structure and sloshing liquid could be made using simplified geometries and the general 
models discussed in chapter 2. 

3.3.2. 1 Model for discrete masses. 

The most straight-forward method of quantitatively determining potential dynamic 
effects of internal structures on the whipping response of a ship is to model the internal 
structure with only a few degrees of freedom (i.e. modes of response), and select 
properties of the structure (i.e. mass, stiffness, damping) which are realistic and will make 
the internal structure respond near resonance (response modes not near resonance would 
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contribute little to the overall response) This amounts to finding a maximum effect 
provided by an internal structure if it were "tuned" to the resonant whipping frequency (or 
the exciting bubble pulse frequency). By using this approach, a quantitative measure of 
the damping effects of internal structures could be obtained for various combinations of 
mass sizes and locations While mass sizes should represent reasonable ship internal 
structures, the locations of the masses may influence any of the low frequency flexural 
modes of the hull, and thus a range of potential effects can be predicted. Figure 3-6 
illustrates a notional internal structure connected to the hull through some combination of 
foundation and mounts with known mechanical properties 




Figure 3-6 Notional internal structure connected to the hull 

The dynamic vibration absorber. For this analysis, a dynamic vibration absorber 
will be used to model a discrete internal structure. The absorber is considered as a single 
discrete mass, as discussed in section 2.3.1, and connected to the ship hull at one of the 
"nodes" defined in section 3.3.1. Ideally, it would be desirable to permit the absorber to 
respond in all six degrees of freedom (three translational and three rotational) However, 
in order to reasonably limit the computational and modeling effort, the absorber will be 
limited to respond only in the vertical or horizontal directions, as shown in figure 3-7. 

This limitation simplifies the computational algorithm, without unduly restricting the value 
of the model. The absorber model can be thought of as being representative of a resonant 
mode of any discrete internal structure. 

Using the approach developed in section 2.3.1 for a two degree of freedom 
dynamic absorber, the equations of motion for the ship hull, as modified for the forces 
exerted on the hull by the absorber may be written in matrix form: 

[MX + CX + KX] iu|] = F,^ + F >lrab « (3-28) 
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where terms on the left hand side are defined by the finite element model for the ship hull 
as described in section 3.3.1. Displacements and forces are vectors positioned at the hull 
nodes, and hull mass, damping, and stiffness are matrices, constructed for the finite 
element formulation. The vector of forces exerted on the hull structure (i.e. at the nodes) 
by the absorber ( F absorber ) is written (ignoring damping associated with the absorber) 

Absorber = { k ( x - X )} = {- mx } (3-29) 

where X (capital) refers to the hull displacement at the node where the absorber acts, and 
x (small case) is the absorber displacement in the direction corresponding to X, as shown 
in figure 3-7. Thus, for a finite element representation of a ship hull with N degrees of 
freedom (i.e. N equations), one additional degree of freedom (i.e. one equation) is 
required for each absorber used 




Figure 3-7. Dynamic vibration absorber free to respond in 
vertical or horizontal directions. 

Suppression of dynamic response of the whipping hull. As was discussed in 
section 2.3.1, the amplitude of motion of a structure could be decreased if the properties 
of a dynamic vibration absorber were "tuned" to match the resonant response frequency of 
the main structure. If it is desired to consider decrease in amplitude of one of the flexural 
modes of hull whipping by application of a dynamic vibration absorber (the lowest flexural 
mode for example), then the hull portion of the two degree of freedom system depicted in 
figure 3-7 could be represented by the modal mass, damping and stiffness for that 
whipping mode (defined by equations (2-13)). 
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In order to quantify the potential dynamic effects of internal structures on the low 
frequency flexural (whipping) modes of a submarine hull, the following process could be 
implemented to develop the equations of motion for the forces exerted by the absorber 
mass (internal structure) on the hull: 

1 Calculate the mode shapes and natural flexural frequencies of hull whipping 
(O n ) by solving the finite element eigenproblem as discussed in section 
3 3 1. These represent the dry' mode shapes and natural frequencies of 
hull whipping. For a submerged submarine, the mass matrix (M) should be 
modified by adding approximate fluid added mass (MJ prior to calculating 
the (modal) natural frequencies since the added mass could have a 
significant effect on the natural frequencies 

2. Calculate the modal mass, damping, and stiffness for the hull whipping 
mode of interest as discussed in section 2.2.2 (equations (2-13)). 

3 Select a "reasonable" absorber mass or mass of an internal structure (m) 
(this could be done for a range of masses). 

4. Calculate the "optimum" absorber stiffness (k), which will "tune" the 
absorber to the resonant frequency of the hull (or the exciting bubble pulse 
frequency): k = m • Q' 

5. Solve for the dynamic response of the structure with the internal absorber 
under specified loading or initial conditions using equations (3-28) and 
(3-29). 

This process will be incorporated into the computational algorithm and used in 
analyses in the next chapter. 

3.3.2. 2 Model of liquid sloshing in a rectangular tank. 

Similarly to the case of internal structure, a quantitative measure of the potential 
damping effects provided by internal sloshing liquids may be made by considering the 
forces exerted on the whipping hull by the modal equivalent mechanical system for the 
complex internal sloshing liquid. In effect, only those sloshing modes likely to be excited 
by the whipping of the hull need be considered in the mechanical equivalent system. The 
sloshing modes not excited may be considered as rigid masses, and thus may be considered 
as part of the main hull mass. 

Using the relations developed in section 2.3.2, the equations of motion for the ship 
hull, as modified for the forces exerted on the hull by a sloshing liquid may be written in 
matrix form: 
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[MX + CX - KX] io( = F {5Mmst + F liquld (3-30) 

where terms on the left hand side are defined by the finite element model for the hull as 
described in section 3.3.1. Displacements and forces are vectors positioned at the hull 
nodes, and hull mass, damping, and stiffness are matrices, constructed for the finite 
element formulation. The vector of forces exerted on the hull structure (i.e at the nodes) 
by the sloshing liquid (F hquid ) can be written (ignoring damping in the tank): 

x x 

F l,qu,d ="Z m n* n = ~ ™ 0 X + X [ K K ~ X)] (3-31) 

n=0 n=l 

where small case letters refer to the mechanical equivalent dynamic parameters for the 
liquid, and subscript n refers to the mode of sloshing. Additionally, the portion of the 
liquid which acts as a rigid body (m 0 ) may be found from equation (2-46) 

x 

m 0 =M l~Z m n ( 3 ’ 32 ) 

n-1 

where M, is the total mass of liquid in the tank. 

These relations may be simplified if it is assumed that only the fundamental 
sloshing mode (n = 1) is excited near its resonance (the higher sloshing modes therefore 
may be considered to act as part of the rigid body mode, m„). Equation (3-31) and (3-32) 
may be rewritten and combined: 

F i,qu,d = -Xm n x n = -(M,X- m,X + m,x,) = (m, - M,)X + [k,(x, - X)] 

n=0 

(3-33) 

Thus, for a finite element representation of a ship hull with N degrees of freedom (i.e. N 
equations), one additional degree of freedom (i.e. one equation) is added for each sloshing 
mode of interest and equations (3-30) and (3-33) form a set of N+l simultaneous 
equations. 

Lateral sloshing in a rectangular tank. Figure 3-8 illustrates a rectangular tank, 
containing a liquid with a free surface, which is oscillating laterally. The potential flow 
solution for the lateral sloshing natural frequencies and hydrodynamic forces are given by 
Vandiver and Mitome 83 for the case where horizontal displacement of the tank is assumed 
to be harmonic and of the form: 

X(t) = Ae‘“ (3-34) 



83 Vandiver and Mitome. op. at.. pp. 26-28. 
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where A is the amplitude of the harmonic motion of the tank wall. The natural frequencies 
are given for mode n: 



K 



co ” = g(2n - 1) — tanh 



2n - 



;rh 



(3-35) 



a v a j 

where g is the acceleration due to gravity, a is the width of the rectangular tank, and h is 
the mean depth of liquid in the tank The modal masses are given for each mode n 



m n = M 



8 tanh 



(2n - 1 )ti 




7i 3 — (2n - 1) 



(3-36) 



a 

As discussed in section 2.3.2, the modal stiffnesses can be computed from the modal mass 
and natural frequency by the expression: 



k„ = co’rn 



(3-37) 




Figure 3-8. Laterally oscillating rectangular tank with standing waves. 

Suppression of dynamic response of the whipping hull. Similarly to the 
application of the dynamic vibration absorber, the motion of the whipping hull could be 
analyzed by considering any onboard tanks as "liquid slosh" dampers. The main 
differences between sloshing liquids and dynamic vibration absorbers discussed in the 
previous section come about because of the numerous sloshing modes which are not 
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excited by the whipping (i.e. the portion of the liquid in the tank which responds as a rigid 
body with the hull, m 0 ). Each sloshing mode which responds may be thought of as an 
individual vibration absorber, but the sloshing modes which do not respond might simply 
be added to the mass of the hull structure. 

Thus, a quantitative feel for the effects of sloshing liquids might be obtained by 
considering the effects of internal structures or absorbers, but also considering the portion 
of the liquids which are not excited, but respond as rigid masses. A more thorough 
discussion of this will be presented in the next chapter. 

3.3.3 External (fluid) forces. 

3.3.3.1 Representation of the hull surface. 

As discussed in section 2.3.1, slender bodies are often represented 
hydrodynamically using a "strip theory" approach since flow around each strip of the hull 
tends to be represented well by considering the flow to be two-dimensional (particularly at 
higher frequencies). Slender submarine hulls, whose length-to-diameter ratios typically 
range from eight to twelve, would therefore be represented well using the hydrodynamic 
"strip theory" approach. 

For purposes of this analysis, a submarine hull will be represented 
hydrodynamically using a number of surface elements ("strips") as shown in figure 
3-9. The centers of pressure of the surface elements will be considered to be located 
coincident with the beam finite element nodes defined in section 3.3. 1.1. Thus, the net 
forces exerted on each of the hull sections by the surrounding fluid could be considered to 
act at the nodes of the finite element model. The properties of each surface element (i.e. 
principle dimensions) would be considered constant for the length of that surface element. 

Because two-dimensional flow is assumed around each surface element, the forces 
on each surface element (and therefore on each finite element node) will generally be 
resolved into two Cartesian coordinate components corresponding to the horizontal and 
vertical directions (normal to the ship axis). Force components on each surface element 
will be computed and applied to the beam finite element node. In the case where a 
particular surface element has appendages, additional contributions to the horizontal and 
vertical forces will be computed for the given appendage geometry. 
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Surface element ('strip') 



Figure 3-9 Hydrodynamic representation of the submarine hull by 
surface elements ("strips") 

3.3.3. 2 Morison relative velocity model. 

As discussed in section 2.3 2, when two dimensional, generally oscillating flow can 
be assumed, and when both fluid inertia and drag forces are important (and wave radiation 
may be ignored), then a Morison's equation may be used to predict hydrodynamic forces 
on a body. It will be assumed for this analysis, that the whipping submarine is sufficiently 
below the free surface such that wave radiation may be ignored. Under this assumption, 
the hydrodynamic forces on each of the surface elements may be predicted using the 
Morison relative velocity formulation given by equation (2-62) and corrected by 
multiplying by the length of the element: 

F hvdro =fp A C m (u-x) + pAx + ^pDC d (u-x)|u-x|jL (3-38) 

where F hydro is the total hydrodynamic force on the surface element, p is the fluid density, 

A is the cross-sectional area of the surface element (normal to the direction of flow), D is 
the characteristic diameter of the surface element (normal to the direction of flow), L is 
the length of the surface element, u is the fluid velocity at the location of the center of 
pressure of the surface element (i.e. the "free field" fluid velocity caused by the pulsating 
explosion bubble and considered relatively unaffected by the presence of the hull), and x is 
the displacement of the center of pressure of the surface element (i.e. the finite element 
nodal displacement) C m and C d are the relative motion fluid inertia and drag 
hydrodynamic coefficients (respectively) and are defined by equation (2-63). 

Determination of hydrodynamic coefficients. For the case of the whipping of a 
submerged submarine, particular care must be taken to select hydrodynamic coefficients 
for analysis which properly represent both the geometry and the conditions of the relative 
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fluid flow around the surface element Generally, it would be expected that separation of 
flow around the hull would occur during submarine whipping motions, particularly around 
those sections with small characteristic dimensions (or containing components with small 
characteristic dimensions). Determination of the extent to which separation occurs and 
the resulting forces generally requires a complex evaluation of the separation and vortex 
shedding mechanisms In the "discrete-point" vortex method 84 , for example, forces on an 
isolated sharp edge may be predicted my modeling the separating shear layers which feed 
the growing vortices. The hydrodynamic coefficients are calculated iteratively over a 
number of cycles by a Fourier averaging process. This type of evaluation, while perhaps 
required for completely rigorous solution of the hydrodynamic forces, is beyond the scope 
of this analysis and will not be considered here 

An initial perspective of the flow fields which could be expected around a 
whipping submarine, and thus of the criteria necessary for initial selection of 
hydrodynamic coefficients for analysis, could be attained by considering past experience 
with ship whipping, and considering typical submarine geometries. As pointed out by 
Hicks 85 , the fundamental (i.e lowest) natural whipping period of most ships is typically 
between 0.3 and 1.0 seconds (fairly well matched with typical bubble pulse periods) For 
typical characteristic dimensions (D) of submarine hull sections (ranges might be between 
0 1 m for small control surfaces, to 15 meters for an entire hull section), characteristic 
frequencies (P) would range from 1 x I0 4 to 7.5 x 10* (higher whipping modes would 
yield even higher (3). These characteristic frequencies correspond to very large Reynolds' 
numbers (Rn) over a large range of Keulegan-Carpenter numbers (KC). 

Extensive experimental studies have been conducted in the marine industry in 
order to obtain values for hydrodynamic coefficients C m , C d , C A , and Cj,, particularly for 
smooth horizontal cylinders, as discussed in section 2.3.2. Even with the extensive 
studies, a great deal of experimental "scatter" exists in the data base, due to such factors 
as free surface effects, surface roughness, and even the random nature and directionality of 
the ocean environment (for field tests). For these reasons, a great deal of care must be 
taken prior to applying coefficient values from the literature to a design or analysis case. 
Additionally, most experimental studies have been carried out at relatively low frequencies 
(and therefore low Reynolds' numbers), and thus there is limited data which could be 
directly applied to submarine whipping without some extrapolation. 



84 Graham. J.M.R.. "The forces on sharp-edged cylinders in oscillatory flow at low Keulegan-Carpenter 
numbers". Journal of Fluid Mechanics. 97. 1980. pp. 331-346. 

85 Hicks. op. cit., p. 393. 
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Extensive laboratory research conducted by Sarpkaya 86 showed that, for smooth 
circular cylinders at high Reynolds' numbers (and therefore high (3), the value of C d 
approaches 0.65 and C m approaches 1.8 (relatively constant over a wide range of KC) 
These coefficients should provide good first approximations for inertia and drag 
coefficients of the circular hull sections (surface elements) for high-frequency whipping. 

Generally, there may be several sharp-edged portions of a submarine surface 
(control surfaces and sail) which may induce significant vortex shedding and therefore 
significant viscous drag on particular hull surface elements Some accounting must 
therefore be made for the additional forces exerted on these surface elements Test data 
on sharp-edged cylinders in high Reynolds' number oscillatory flow is generally lacking in 
the literature. Some data is presented by Bearman. et. al. 87 for flat plates in oscillatory 
flow at Keulegan-Carpenter numbers between 1 and 10 (relatively low) which might be 
extrapolated to predict that C d would fall in the range of 4-6 and C m would fall in the 
range of around 1.0 for higher Reynolds' number (higher frequency) flows. While these 
values are probably not experimentally sound, they should at least provide reasonable first 
order approximations of additional viscous forces exerted on certain hull surface elements 
by sharp-edged control surfaces or sail 

The modeling of hydrodynamic forces on a whipping submarine using equation 
(3-38) will be included as part of the computational algorithm developed in section 3 4 
One analysis application will be made in chapter 4 where whipping response of the 
submarine is due only to an initial displacements (i.e. free whipping), and therefore no 
"free field" fluid velocity or acceleration terms (i.e. u and u) will be included, and the 
hydrodynamic forces will be due only to the motion of the hull in otherwise still water 
Additional analyses will include cases where the submarine is excited by a pulsating 
explosion bubble, and "free field" fluid motion will be included in the hydrodynamic 
model. 

3.3.3.3 Model of fluid loading from a pulsating explosion bubble. 

In analysis cases where it is desirable to couple ship whipping with the loading 
imposed by the pulsation of an explosion bubble, it is necessary to define a time-varying 
formulation for the "free field" fluid velocity and acceleration associated with the pulsating 



86 Sarpkaya. T.. "In-line and transverse forces on cylinders in oscillating flow at high Reynolds' numbers". 
Proceedings of the Eight Offshore Technology Conference. OTC 2533. 1976. pp. 95-108 
87 Bearman. P W.. Dow me. M.J.. Graham. J.M.R . and Obasaju. ED. "Forces on cylinders in viscous 
oscillator, flow at low Keulegan-Carpenter numbers". Journal of Fluid Mechanics. 154. 1985. pp. 344- 
345. 
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bubble (i.e. u and u) Because derivation of equations for bubble dynamics is a lengthy 
process and not a goal of this analysis, formulation available in the literature will be only 
briefly summarized A particularly useful formulation for spherically symmetric, pulsating 
and migrating explosion bubbles is presented by Wilkerson 88 , which is based upon the 
published theories of Herring 89 , Friedman 90 , Taylor 91 , Cole 92 , and Hicks 93 This 
formulation includes accounting of effects due to the presence of a free surface and drag 
forces on the migration of the bubble. 

The free field equations of motion for a spherically symmetric, pulsating and 
migrating explosion bubble are described adequately using inviscid. irrotational fluid flow 
theory (i.e. potential flow). Figure 3-10 illustrates an appropriate coordinate system which 
encompasses the explosion bubble, its "image bubble" (to account for the presence of the 
free surface), and an arbitrary point in the fluid. The fluid velocity potential function (<j>) 
resulting from a pulsating and migrating explosion bubble may be written as the 
superposition of potentials resulting from a simple source and its image (representing 
pulsation of the spherical bubble) and potentials resulting from a dipole and its image 
(representing migration of the spherical bubble) 

4) = — + cos0, - — + ^|-cos0, (3-39) 

r i r f r : r : 

where e, and e ; are the time-dependent strengths of the sources and dipoles 
(respectively), r, , r ; are the radial distances from the center of the bubbles to the arbitrary 
point in the fluid, and 0, ,0, are the angles from the vertical (as shown in figure 3-10). 

The first term on the right hand side is the potential due to the source representing the 
pulsation of the explosion bubble. The second term is the potential due to the dipole 
representing the migration of the explosion bubble. The last two terms represent the 
potentials due to the effects of the "image bubble" which accounts for the effect of the free 
surface (i.e. satisfies the boundary condition at the free surface). 



88 Wilkerson. S. A.. "Elastic Whipping Response of Ships to an Underwater Explosion Loading". Master of 
Science in Engineering Thesis. George Washington University. 1985. 

89 Hemng. C.. "Theory of the Pulsations of the Gas Bubble Produced by an Underwater Explosion", 
Underwater Explosions Research, I'ol. II - The Gas Globe. Office of Naval Research, 1950. 

90 Friedman. B.. "Theory of the Underwater Explosion Bubbles". Underwater Explosions Research, I 'ol. II 
- The Gas Globe. Office of Naval Research. 1950. 

9l Taylor. G.I.. "Vertical Motion of a Spherical Bubble and the Pressure Surrounding It", Underwater 
Explosions Research, I'ol. II - The Gas Globe. Office of Naval Research. 1950. 

92 Colc. R.H.. Underwater Explosions. Princeton University Press. Pnnceton. NJ. 1948. 

93 Hicks. A.N.. "The Theory of Explosion-Induced Ship Whipping". Xaval Construction Research 
Establishment, XCRE Report R579. 1972. 
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"Image" bubble 




Figure 3-10. Coordinate system for explosion bubble dynamics. 
Explosion and "image" bubble representation. 



The time-dependent source and dipole strengths are those required to maintain 
continuity of the normal fluid particle velocity at the bubble surface (i.e. the normal fluid 
particle velocity must equal the rate of change of the bubble radius, a , corrected for the 
migration velocity). The required source and dipole strengths are: 

e, = a : a 



.3 



e, - 



v m -■ 



a'a 

4c! 7 



(3-40) 



where v m is the upward velocity of the (migrating) bubble and d is the depth of the 
bubble. Using conservation of energy, the differential equations which govern the 
pulsation and migration of the bubble (i.e. solve for the bubble radius and location) may be 
derived 94 . The differential equations (in non-dimensionalized form) which describe bubble 



94 Wilkerson. op. cit.. pp. 7-24. 



77 



pulsation and migration are presented as pan of figure 3-11, along with other defined 
constants and non-dimensionalized variables used in the solution of the bubble dynamics. 

The first order differential equations which govern the pulsation and migration of 
the bubble may be integrated using a time-stepping integration procedure such as a Runge- 
Kutta method It is recommended 95 that a variable time-stepping procedure be used 
because of the rapid changes in radius during the phases of bubble minima. With the 
calculated bubble radii and positions (with time), the velocity potential may be calculated 
using equations (3-39) and (3-40), and the "free field" fluid velocity for any arbitrary point 
in the fluid may be calculated by differentiating the velocity potential (equation (2-47)) 
Wilkerson 96 derived expressions for fluid velocities and accelerations (horizontal and 
vertical components) at an arbitrary point in the fluid as functions of position and source 
and dipole strengths. The resulting expressions (given in Cartesian coordinates of figure 
3-10) are given in figure 3-12. 



95 Hicks. op. cit.. p 159 
96 Wilkerson. op. cit.. pp. 25-29. 
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Differential equations which govern bubble pulsation and migration: 



a - a 



a = - 




c r 'f 1 2 pa 



a v 35 




( 7 - 1 ) k (3a * a/. 



45 : 
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Initial Conditions: a = k 4 5 1^—7- 



,, = -d / L a = 0 /. = 0 



Non-dimensionalized defined variables: a = a / L (non-dimensional bubble radius) 



D 0 - d + 33 d is initial depth of charge W is charge weight of TNT 
Control variables: s (0 for non-migrating; 1 for migrating) 

(3 (0 for no free-surface effect, 1 for free-surface effect) 

Figure 3-11 Non-dimensionalized variables and differential equations which describe 
the pulsation and migration of an underwater explosion bubble. 



; = y b / L (non-dimensional bubble depth) 
x = t / T (non-dimensional time) 



where: L a 13.6(W/ D 0 )‘ 0 for TNT 



k * (0.0552 )Dq : 5 for TNT 
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y * 1.25 for TNT C D = 2.25 
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where source and dipole strengths may be defined (in non-dimensionalized form) 

L J a 3 A 
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and: v m = -^L/T = -AL/T (the upward velocity of the bubble) 



Figure 3-12. Fluid velocities and accelerations (in relative Cartesian coordinates) for 
given source and dipole strengths and upward bubble velocity 



80 



3.4 Formulation of the computational algorithm. 

In order to integrate each of the desired force components presented in the 
previous sections, a comprehensive computational algorithm must be developed. The 
required algorithm must be capable of performing the following: 

1. Compute the time histories of the explosion bubble radius and depth by 
integration of the differential equations governing bubble dynamics using a Runge- 
Kutta time integration procedure. 

2. Compute the "free-field" fluid velocities and accelerations due to the 
pulsating and migrating explosion bubble These velocities and accelerations must 
be calculated at the locations of each of the finite element nodes used to model the 
dynamics of the hull structure at each iteration at each tune step , and must be 
resolved into the vertical and horizontal directions defined by the ship coordinate 
system (a change of coordinates from the bubble coordinate system to the ship 
coordinate system is required). 

3. Compute the external and internal loading on the ship hull using the 
appropriate models developed in the previous sections. 

4. Iteratively determine the hull structure displacements, velocities, and 
accelerations using the Newmark method (with Newton-Raphson iterations 
conducted at each time step) as described in section 3.2. Also, iteratively 
determine displacements, velocities, and accelerations of internal masses as 
required. 

A computational algorithm satisfying the above requirements may be generated 
using a suitable computer language, separating the algorithm into sections or subroutines 
appropriate for computation using the chosen language. For this application, the 
algorithm will be implemented on a personal computer using subroutines written for the 
MATLAB® computer program Appendix A presents MATLAB® subroutines written 
for this analysis, each of which addresses specific portions of the overall algorithm. The 
subroutines are presented as follows: 

a. Subroutine IN.M - Serves as an input file for other subroutines. Data is loaded 
from formatted batch-type files. Arrays and constants used in other subroutines are 
generated 

b Subroutine HULL M - Calculates global mass and stiffness matrices for the hull 
structure, solves the generalized eigenproblem to obtain the "dry mode" natural 
frequencies and mode shapes, calculates an appropriate viscous hull damping matrix for 
the hull using a Rayleigh damping approach, and calculates vectors for modal mass, 
damping, and stiffness. 
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c. Subroutine WETMODE. M - Calculates an approximate fluid added mass 
matrix for the hull using a Morison equation approach (i.e. strip theory), adds the fluid 
added mass matrix to the dry hull mass matrix and solves the generalized eigenproblem to 
obtain the "wet mode" natural frequencies and mode shapes, and calculates vectors for 
"wet mode" modal mass, damping and stiffness. 

d Subroutine FREEWHIP M - Calculates nodal point response (displacements, 
velocities, and accelerations) during free vibration with initial displacements given by 
scaled mode shape vectors for specified modes Response is calculated using the 
Newmark time integration method with iteration conducted at each time step using a 
Newton-Raphson method. This subroutine is designed primarily to compare damping 
associated with hull mechanisms (i.e. material hysteretic and structural) and damping 
associated with external hydrodynamic forces 

e. Subroutine BUBBLE. M - Computes and saves the time histories of the 
explosion bubble radius and bubble depth by integration of the differential equations 
governing bubble dynamics. Time histories are computed using a Runge-Kutta integration 
method with a variable time step. The differential equations governing bubble radius are 
called from function routines RADFN.M and RFUNC.M The differential equations 
governing bubble depth are called from function routines DPTHFN.M and DFUNC.M. 

f. Subroutine BUBWHIP.M - Calculates the nodal point response (displacements, 
velocities, and accelerations) during forced vibration due to a pulsating explosion bubble. 
Response is calculated using the Newmark time integration method with iteration 
conducted at each time step using a Newton-Raphson method. Additionally, local axial 
strains are calculated at specified hull locations using simple beam theory (used for 
comparison of the computational algorithm with model test data) The free-field fluid 
velocities and accelerations due to the pulsating explosion bubble are calculated at each 
iteration when function routine FLUID. M is called from within BUBWHIP.M. This 
subroutine is designed primarily to provide a comprehensive computation of explosion- 
induced whipping effects Additionally, it may be used to quantify damping associated 
with internal structures and to provide for comparison of the computational algorithm with 
model test data. 
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4. ANALYSIS OF WHIPPING OF A SUBMERGED SUBMARINE 



4.1 Introduction. 

In this chapter, the prediction of the whipping response of a submerged submarine 
is discussed, including an analysis of specific damping mechanisms The computational 
algorithm developed in the previous chapter is utilized to predict the response of a 
notional 7400 long ton fast attack submarine (SSN). Additionally, the computational 
algorithm is used to predict the response of the U S Navy test platform "Red Snapper", 
subjected to specified explosive charge weight/standoff combinations for comparison with 
model test data. The latter analysis is presented in non-specitic response units due to 
security concerns with the actual test data used for comparison. 

General characteristics of the notional submarine (SSN) are given in Table 4-1. 
From these general characteristics, specific characteristics (beam element data, hull section 
data, etc.) are formulated and provided as input via subroutine IN.M. The "dry" mode 
natural frequencies and mode shapes are calculated using subroutine HULL.M, and the . 
"wet" mode natural frequencies and mode shapes are calculated using subroutine 
YVETMODE. M (i.e. the effects of fluid added mass are included) Table 4-2 gives both 
"dry mode" and "wet mode" natural frequencies computed for the first 12 displacement 
modes. It should be noted that this represents the first 6 displacement modes in each 
plane (i.e. horizontal and vertical), of which the first 2 are considered rigid hod}' modes 
(corresponding to zero natural frequencies), and the remaining 4 modes are considered as 
the fundamental flexural modes. The horizontal "wet" mode shapes are plotted (as 
normalized nodal displacements) and shown in figure 4-1 (rigid body mode shapes) and 
figure 4-2 (fundamental flexural mode shapes). 



Table 4-1. General characteristics of the notional SSN 



Length (ft) 


350 


Displacement/weight (long tons) 


7440 


Maximum beam/diameter (ft) 


35 


Hull material 


HY-80 steel 


Number of hull sections/beam elements 


20/19 


Appendages 


Sail, bow planes, stern planes, rudder 



83 



Table 4-2. Rigid body and fundamental flexural "dry" mode and "wet" mode natural 
frequencies for notional SSN. 



Mode 


"Dry mode" natural frequency, © 
(hz, rad/sec) 


"Wet mode" natural frequency, co 
(hz, rad/sec) 


1 Rigid body 


0 


0 


2 Rigid body 


0 


0 


3 Rigid body 


0 


0 


4 Rigid body 


0 


0 


5 Flexural (V) 


1.984, 12.466 


1.370, 8 608 


6 Flexural (H) 


1.984, 12.466 


1.373, 8 627 


7 Flexural (V) 


4.498, 28.262 


3.060, 19.227 


8 Flexural (H) 


4.498, 28.262 


3.074, 19.315 


9 Flexural (V) 


7.054, 44.322 


4.827, 30.329 


lOFlexural (H) 


7.054, 44.322 


4.863, 30.555 


1 1 Flexural (V) 


9.802, 61.588 


6.907, 43.398 


12Flexural (H) 


9.802, 61.588 


6.922, 43.492 



Wet mode shape of rigid body mode in horizontal (y=0) plane 
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Wet mode shape of flexural mode 1 in horizontal (y=0) plane 




Wet mode shape of flexural mode 2 in horizontal (y=0) plane 




Wet mode shape of flexural mode 3 in horizontal (y=0) plane 




Wet mode shape of flexural mode 4 in horizontal (y=0) plane 




Figure 4-2. Horizontal fundamental flexural "wet" mode shapes for notional SSN. 
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4.2 Prediction of "dry" whipping response. Hull damping. 

The "dry" whipping response of the SSN is computed for given initial conditions 
using subroutine FREEYVHIP M (with fluid density set to zero in Morison’s equation). 

For free "dry" whipping (as well as free "wet" whipping discussed subsequently), initial 
nodal point displacements are specified as scaled horizontal mode shape vectors, and time 
histories of nodal point displacements, velocities, and accelerations are computed based 
upon these initial displacements. Thus, initial loading on the hull structure is due onlv to 
internal forces due to structural stiffness. 

For this analysis, initial displacements are set to the first and second horizontal 
fundamental flexural mode shapes (scaled by 2% of the total ship length). Appendix B 
presents plots of horizontal displacements, velocities, and accelerations for nodes 7, 10, 
and 20 for initial displacements (a) scaled from the 1st horizontal fundamental flexural 
mode shape (i.e. mode 6 in Table 4-2), and (b) scaled from the 2nd horizontal fundamental 
flexural mode shape (i.e. mode 8 in Table 4-2). 

As discussed in chapter 3, it has been assumed that the inherent hull damping (i.e. 
due to material hysteresis and structural damping) would reasonably follow experimental 
data for steel cylindrical marine structures and ship masts Thus, it was assumed that 
modal damping ratios on the order of 1% for hull damping would be sufficient and 
Rayleigh coefficients could be calculated on this basis Thus, reductions in free response 
amplitude with time (damping) seen in graphs presented in Appendix B are merely 
manifestation of this initial assumption, and no further conclusions may be drawn 
concerning the validity of this initial assumption. 

However, an approximate graphical check could be performed to verify that 
damping in the case of "dry" whipping is indeed represented computationally by the 
assumed 1% modal damping. Indeed, measurement of amplitude log decrements (recall 
equation (1-1)) and calculation of modal damping ratios (recall equation (2-40)) for both 
mode 1 horizontal flexure (represented by node 10 displacement) and mode 2 horizontal 
flexure (represented by node 7 displacement), show that graphically determined modal 
damping ratios are approximately 1%. 

4.3 Prediction of "wet” whipping response. Hydrodynamic damping. 

The "wet" whipping response of the SSN is computed for the same given initial 
conditions as the "dry" mode analysis, also using subroutine FREEWHIP.M (with 
inclusion of hydrodynamic forces by Morison's equation). Appendix C presents plots of 
horizontal displacements, velocities, and accelerations for nodes 7, 10, and 20 for initial 
displacements (a) scaled from the 1st horizontal fundamental flexural mode shape, and (b) 
scaled from the 2nd horizontal fundamental flexural mode shape. 
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Figure 4-3 shows comparisons between nodal point displacements for "wet" and 
"dry" cases for initial displacements scaled from mode 1 and mode 2 flexural mode shapes 
Even though modal damping cannot be directly applied to the nonlinear hydrodynamic 
damping forces, a calculated "equivalent" modal damping might still serve as a measure of 
hydrodynamic damping in comparison to the "dry" case As was done for "dry" whipping, 
a graphical approximation of "equivalent" modal damping ratios can be made for mode 1 
horizontal flexure and mode 2 horizontal flexure (represented by node 10 and node 7 
displacements, respectively). As expected, the inclusion of external hydrodynamic forces 
(represented by Morison's equation) provide for increased modal damping ratios for both 
flexural modes. Both mode 1 and mode 2 flexure are characterized by "equivalent" modal 
damping ratios of approximately 2.2%. 



Horizontal displacement, node 10 





Figure 4-3. Comparison between nodal point displacements for "dry" and "wet" hull 
flexure. Top: mode 1 flexure, represented by node 10 displacement. 
Bottom: mode 2 flexure, represented by node 7 displacement. 
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4.4 Prediction of whipping response due to explosion bubble pulse loading. 

The previous sections presented analysis based upon an assumption of initial 
displacements set to flexural mode shapes (i.e. free whipping), aimed only at providing 
some quantitative measure of hydrodynamic damping. As discussed in chapter 1, the main 
concern with submarine whipping involves the quasi-periodic bubble pulse loading 
resulting from near to intermediate stand-off explosions It is this quasi-periodic bubble 
pulse loading which can reinforce the free vibration and magnify the whipping response, 
particularly when the period of the bubble pulse is in the vicinity of one of the natural 
periods of hull flexure Indeed, if the overall dissipative forces (damping forces) are low, 
whipping response magnitudes (and therefore resulting hull damage) could be quite 
substantial. 

To see potential dynamic effects associated with matched bubble pulse (i.e. bubble 
pulse frequency matched to ship's flexural frequency), explosive charge characteristics 
(charge weight, depth, etc.) must first be defined Using the formulation for bubble 
pulsation and migration discussed in chapter 3, and the known "wet" natural frequencies 
of hull flexure (Table 4-2), explosive charge characteristics could be iteratively 
determined. Table 4-3 specifies explosive charge characteristics used in this analysis to 
match bubble pulse frequency with 1st and 2nd horizontal flexural mode natural 
frequencies (i.e. modes 6 and 8 in Table 4-2). A charge depth of 100 feet was selected, as 
deeper depths require unrealistically large charges to match the 1st fundamental flexural 
frequency (1.37 Hz). Ship location relative to the charge was chosen to promote 
excitation of the mode (both frequency and shape) of interest. 

Figure 4-4 shows plots of bubble radius, depth, and migration velocity for the 
explosive charge weight/depth combination selected to excite mode 2 horizontal whipping 
of the SSN. It can be seen from the plot of radius vs. time that the bubble period is not 
constant with time, but increases as the bubble migrates toward the surface This is an 
important phenomenon, which limits the ability of the pulsating bubble to excite a 
particular flexural natural frequency for a prolonged period of time (thus, the bubble pulse 
could only be considered to be quasi-periodic in nature) 

The whipping response of the SSN is predicted for the cases when excited by 
explosion bubble pulse characterized by charges specified in Table 4-3. Subroutine 
BUBBLE. M predicts explosion bubble radius and depth time histories for each of the 
charge weight/depth combinations selected Subroutine BUB WHIP. M predicts the nodal 
point response of the SSN, initially at rest and subject to the incident bubble pulse loading. 
Appendix D presents plots of displacements, velocities, and accelerations for nodes 7, 10, 
and 20 for the charges designed to excite horizontal modes 1 and 2 flexural whipping 
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Table 4-3. Charge characteristics selected to excite mode 1 and mode 2 horizontal 
flexural whipping for SSN. 





First horizontal 
flexural mode 


Second horizontal 
flexural mode 


Charge weight (lb TNT) 


1000 


85 


Charge depth (ft) 


100 


100 


Bubble pulse frequency (Hz) 


1.35 


3.05 


Depth of ship (ft) 


70 


80 


Charge location along length of 
ship (ft) 


157.5 (node 10) 


105.0 (node 7) 


Charge standoff from ship 
centerline (ft) 


33.5 


15 


Charge depth below ship 
centerline (ft) 


30 


20 



Because of the swiftness with which the (large) bubble migrates to the surface 
from its initial depth for the larger/lower frequency charge (mode 1 whipping), the 
calculated response time period was limited to the first 2.5 seconds. Despite the short 
time of excitation (only 3 flexure/bubble periods), it can be easily observed that the quasi- 
periodic bubble pulse does indeed magnify the whipping response For the smaller/higher 
frequency charge (mode 2 whipping), the magnification effect can be more easily observed 
over a larger number of whipping cycles. 
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Bubble radius vs. time 




Bubble depth vs. time 




Time (sec) 



Upward bubble velocity vs. time 




Figure 4-4. Plots of bubble radius, depth, and migration velocity for the explosive charge 
weight/depth combination selected to excite mode 2 horizontal whipping of the SSN. 
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4.5 Damping from discrete internal structures and sloshing liquids. 

Discrete internal structures. As discussed in chapters 2 and 3, whipping energy 
can be transferred to internal structures within the hull, particularly when the internal 
structure is "tuned" to the bubble pulse frequency (and one of the natural frequencies of 
hull flexure). Using the model developed in chapter 3 (which models the internal structure 
as a discrete mass or dynamic vibration absorber), a brief analysis is conducted to gain 
some insight into potential dynamic effects of these internal structures on the whipping 
response of a hull. 

The analysis discussed in the previous section is extended to include "tuned" 
vibration absorbers (internal masses on springs) where each absorber is "tuned" to the 
matched bubble pulse/w hipping frequency The effects of a "tuned" internal absorber are 
predicted for the case of the bubble pulse matched to the 2nd horizontal flexural 
frequency. The internal absorber is located at node 7 (i.e. "attached" to node 7) and free 
to move in the horizontal direction only, as discussed in section 3.3. Thus, the effect of 
the absorber is on the 2nd horizontal whipping mode 

The effect of a 200 long ton internal mass, "tuned" to the 2nd horizontal whipping 
mode and "attached" to the hull at node 7, is considered. Appendix E presents plots of 
horizontal displacements, velocities, and accelerations for nodes 7, 10, and 20, as well as 
horizontal displacement, velocity, and acceleration of the absorber mass. Figure 4-5 plots 
horizontal displacement of node 7, with and without the 200 ton internal mass. 

Several interesting phenomenon are evident from the plot and should be briefly 
noted. First, the presence of the absorber has a significant overall effect on reducing the 
magnitude of response. While no specific measure of "damping" can be made, there is 
obviously significant energy transferred into the absorber. Second, the ability of the 
absorber to "absorb" whipping energy is very sensitive to the frequency of excitation. As 
discussed in the previous section, the bubble period (and therefore the frequency of 
whipping excitation) is not constant , but rather increases as the bubble migrates toward 
the surface By "tuning" the internal absorber to the 2nd horizontal flexural frequency, the 
absorber is not then truly tuned to the excitation frequency (recall that matching bubble 
period to hull flexural frequency was accomplished by iteration of average bubble periods 
- over the 3 seconds of bubble pulsation). Thus, the usual "internal absorber effect" (i.e. 
reducing vibration amplitude of the main structure to zero) would not occur. Figure 4-6 
plots the horizontal displacement of node 7 and the horizontal displacement of the 
absorber mass. Because the oscillation of node 7 and the oscillation of the absorber mass 
occur at slightly different frequencies, a sharp (but temporary) reduction in node 7 
amplitude can be seen around 2.5 seconds, when the absorber shifts from a lag angle to a 
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lead angle (i e the node and the mass are temporarily 180° out of phase) What can be 
deduced from this is that internal structures do have potential for significant absorption of 
hull whipping energy, even though bubble pulse frequencies are not truly constant. 
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Liquid sloshing. As discussed in chapters 2 and 3, the fundamental mechanism by 
which energy is "transferred" to internal liquids through the mechanism of sloshing can be 
thought of much in the same way as that by which it is transferred to internal structures 
In fact, liquid sloshing may be modeled quite accurately using a mechanical equivalent 
system approach, where each sloshing mode is considered to have its own modal mass, 
damping and stiffness, each responding as its own dynamic system (i.e. its own dynamic 
absorber). Using the formulation for natural sloshing frequency, modal mass, and modal 
stiffness presented in chapter 3 for lateral sloshing in a rectangular tank, the discussion of 
damping through liquid sloshing follows quite nicely behind that of the discrete internal 
structure. Recalling equations (3-35) through (3-37), it is clear that the natural frequency, 
modal mass, and modal stiffness of each mode of sloshing are dependent only upon the 
mode number (n), the total mass of liquid in the tank (M,), and the geometry of the tank 
(i.e. tank depth and width) 

Using the standard rectangular tank formulations, it is a trivial procedure to 
generate plots of natural frequency vs. sloshing mode and sloshing modal mass vs 
sloshing mode. When this is done over a wide range of possible tank configurations, it 
quickly becomes clear that liquid sloshing could not become a significant contributor to 
the dissipation or absorption of hull whipping energy. Figure 4-7 plots both natural 
sloshing frequency (co n ) and sloshing modal mass (m n ) for a 20' x 20' x 20' rectangular 
tank (a total liquid weight of 228 5 long tons). Two interesting things may be noted for 
this very large tank (much larger than any single, unbaffled tank on any submarine). First, 
the only sloshing modal mass of any substance (i.e. which might effect hull whipping) is 
that of the first/fundamental sloshing mode (i.e. 60 long ton). Second, the natural 
frequency of this fundamental sloshing mode is only 2.24 rad/sec (0.35 Hz), well below 
any natural flexural frequency. As other tanks are considered, it is found that the 
fundamental sloshing modes (of all reasonable tank sizes) fall well below 1 Hz. While 
higher sloshing modes may have natural sloshing frequencies in the range of hull whipping 
frequencies, the modal masses of these modes are significantly smaller than would 
contribute to any noticeable reduction in whipping response. Appendix F presents plots of 
natural sloshing frequency vs. sloshing mode and sloshing modal mass vs. sloshing mode 
for several representative rectangular tank geometries. 

From this basic analysis of natural sloshing frequencies and sloshing modal masses, 
it can be justifiably deduced that liquid sloshing would have little impact on the whipping 
response of most ships. Because the significant sloshing modal masses have natural 
frequencies well below hull flexural frequencies, it can also be deduced that the major 
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effect of liquids with free surface would only be toward the addition of total ship mass (i.e 
the significant sloshing modes would act only as rigid bodies) 
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Figure 4-7. Natural sloshing frequency and sloshing modal mass vs. mode number for 
lateral sloshing of a (20 1 x 20' x 20') rectangular tank. 



4.6 Comparison of computational algorithm with the "Red Snapper" model test. 

The computational algorithm for submarine whipping developed in chapter 3, and 
applied to the whipping of a notional SSN in the previous sections, is used to predict the 
response of the U.S. Navy test platform "Red Snapper". Thus, some "verification" of the 
computational algorithm is made with actual model test data. Due to the security 
concerns with the actual test data used for comparison, the following is presented in non- 
specific response units (i.e. the y-axes are unlabeled). 
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The "Red Snapper" model is on the order of 1/4 the size of the notional SSN 
discussed previously (geometrically), and thus natural frequencies of hull flexure are much 
higher than the notional SSN. Whipping response was measured in terms of local axial 
strains (at specific beam element locations) and horizontal and vertical velocities (at 
specific node locations). For comparative purposes, subroutine BUBWHIP.M includes a 
scheme for calculating local axial strains at specified finite beam elements This scheme is 
based upon simple beam theory an provides local axial strains at port and starboard 
extreme fibers, as well as crown and keel extreme fibers for the specified elements (see 
subroutine BUBWHIP.M in Appendix A). 

Three specific model tests, which were conducted on the "Red Snapper" model, 
are used for comparison One test was designed such that the explosion bubble pulse 
would excite mode 1 horizontal whipping. A second test was designed such that the 
bubble pulse would excite mode 1 whipping in both horizontal and vertical planes A third 
test was designed such that the bubble pulse would excite mode 2 horizontal whipping. 

The design of the bubble pulse to match the flexural frequencies and mode shapes was 
carried out in a manner similar to that discussed in section 4.4 In all cases, there where 
no internal structures or liquids, and thus all forces were either hydrodynamic or structural 
in nature. 

General comparisons. Figure 4-8 shows starboard-side fiber strain histories for 
elements 6 and 9 for Test 1 (horizontal mode 1 whipping), and strains calculated using the 
computational algorithm In general, it appears that calculated strains match well with 
measured through the 2nd bubble period, but then over predict strains. Figure 4-9 shows 
horizontal velocity histories for nodes 10 and 20 for Test 1, and velocities calculated using 
the computational algorithm. Calculated velocities for this test are only slightly under 
predicted. 

Comparisons between measured and calculated strains and velocities for Test 2 
(mode 1 horizontal and vertical whipping) and Test 3 (mode 2 horizontal whipping) are 
presented in Appendix G. In general, the computational algorithm slightly under predicts 
and predicts slightly early for both strains and velocities for Test 2. For Test 3, the 
computational algorithm predicts quite well through the first two bubble periods, but then 
grossly over predicts both strains and velocities 
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Strain Strain 



Axial strain vs. time, element 6 




Axial strain vs. time, element 9 




Time (sec) 



Figure 4-8. Starboard-side fiber strain histories for elements 6 and 9 for Test 1 of 
the "Red Snapper" model test (horizontal mode 1 whipping), and strains calculated using 

the computational algorithm 
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Velocity 



Horizontal velocity vs. time, node 10 




Horizontal velocity vs. time, node 20 




Time (sec) 



Figure 4-9 Horizontal velocity histories for nodes 10 and 20 for Test 1 of 
the "Red Snapper" model test (horizontal mode 1 whipping), and velocities calculated 

using the computational algorithm 
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Discussion of general comparison. There are several possible explanations for the 
failure of the computational algorithm to properly predict strains and velocities for the 
"Red Snapper" tests 

a. The bubble model used for formulation of equations governing bubble dynamics 
has been known to be inaccurate past the 2nd bubble period 97 . Although no satisfactory 
purely physical model exists for long-time bubble dynamics (some more restrictive 
empirical relations are available), this particular bubble model has proven to be fairly 
accurate through the 2nd bubble period The reasons for the inaccuracy past the 2nd 
bubble period lie in the inability of the model to properly account for dissipation of bubble 
energy during bubble pulsation and migration. Thus, the model believes that the bubble 
maintains more of its original energy throughout the longer time history This translates 
into greater bubble radii and greater velocities being modeled in the computational 
algorithm (resulting in over predicting hydrodynamic loading past the 2nd bubble period). 
Thus, some explanation for overproduction of response past the 2nd bubble period exists. 

b. The use of the hydrodynamic "line-structure" model in the computational 
algorithm (i.e. hydrodynamic loading being applied at the nodes along the ship centerline, 
vice on the outer hull of the ship), leads to inaccuracies in hydrodynamic loading. This is 
particularly important when the explosive charge is close aboard In this case, the actual 
hydrodynamic loading is more appropriately seen directly on the outer hull , while the 
computational algorithm believes the loading to be applied at the ship's centerline (i.e. 
calculated "ffee-field" fluid velocities and accelerations). Thus, an error is introduced due 
to the unaccounted-for distance between the centerline and the outer hull of the ship. For 
large ships, subjected to charges with large standoffs, this geometric inaccuracy is less 
important. But, for smaller ships subjected to near-standoff charges, this inaccuracy could 
lead to significant error. 

c. There is a general inability of Morison's equation to adequately account for the 
transient-type vortex shedding and convection (viscous drag forces), which is undoubtedly 
important in this application. The original Morison's equation was designed to account for 
hydrodynamic loading for the relatively low-frequency wave loading on pile structures in 
regular sea waves (i.e. quasi-steady). There is no adequate data in the literature to 
support extension of the Morison formulation to transient-type analyses. This is 
particularly important in the justification and selection of the coefficients used in Morison's 
equation. As discussed in chapter 3, coefficients for the computational algorithm were 
merely extrapolated from low frequency data. Additionally, this data was itself obtained 



97 \lisovec. A.. "Comments on SITE meeting of 5/10/94". (Unpublished meeting minutes). 1994. pp. 1-2. 
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by Fourier time-averaging of measured forces (i.e. quasi-steady). For transient-type 
analyses, such as submarine whipping, a more complex approach would probably be 
appropriate to properly account for the transient vortex shedding and convection. 
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5. CONCLUSIONS 



An analytical approach has been presented which predicts the elastic whipping 
response of a submerged submarine to the pulsating bubble of a nearby underwater 
explosion. The approach presented concentrates on mechanisms associated with energy 
dissipation (damping), and includes specific mechanisms associated with hull damping 
(material and structural), internal damping (internal structures and sloshing liquids), and 
external damping (hydrodynamic damping forces in conjunction with traditional inertial 
loading forces). The formulation utilizes a numerical time-domain solution technique to 
directly solve the nonlinear equations of motion. 

While the general mechanisms associated with damping are reasonably well 
understood, application of this understanding to complex high-frequency, non-harmonic 
vibration problems (such as ship whipping due to explosion bubble pulse loading) has not 
been met with much success. For this reason, the approach has been designed to utilize a 
robust solution technique (the time-domain solution), and incorporate physically-based 
models where possible. 

Modal damping ratios for material and structural damping were assumed to be on 
the order of 1 % of critical for fundamental flexural modes, based upon data available 
from steel marine structures. Thus, results for “dry” whipping only reflect this initial 
assumption. With this initial assumption for material and structural damping, “equivalent” 
modal damping ratios for mode 1 and mode 2 whipping (the first two fundamental flexural 
modes) for the case where hydrodynamic damping effects are included (via Morison’s 
equation), are computed to be on the order of 2.2 %. 

The effects of internal structures and sloshing liquids on whipping response is 
approximated using a 200 long ton “tuned” mass damper, located to effect mode 2 
whipping. The effects, illustrated in figure 4-5, are significant only when bubble pulse 
frequencies are precisely equal to the “tuned” frequency of the internal structure. The 
effects of internal sloshing liquids are shown to be negligible because of the low natural 
fundamental sloshing frequency associated with larger liquid tanks. 

Several of the selected modeling techniques are shown by comparison with actual 
test data to be questionable. Firstly, the representation of hydrodynamic forces solely by a 
Morison's equation approach fails to fully account for the expected vortex shedding and 
convection forces during the transient vibration. The selection procedure (or more 
appropriately, the calculation procedure) of hydrodynamic coefficients appropriate for ship 
whipping must be better understood and implemented. Mere extrapolation of known data 
based upon harmonic, low frequency motion of small scale cylinders and flat plates is 
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probably inaccurate. A more robust model for hydrodynamic forces, coupling numerical 
simulation with experimental verification should be developed. 

A second modeling technique of questionable success is the bubble pulsation and 
migration equations used to produce the free-field fluid velocities and accelerations. 
Because the model does not account for proper dissipation of energy by the pulsating 
explosion bubble, it is believed that the model over predicts fluid loading after the 2nd 
bubble period. A more robust model for fluid loading from the pulsation of the explosion 
bubble, again, coupling numerical simulation with experimental verification should be 
more successful. 

Finally, the general nature of the line-structure model for the ship structural and 
hydrodynamic interaction (i.e. the fluid-structure interaction) introduces geometric errors, 
particularly when the explosion is close aboard. A more robust, 3-dimensional fluid- 
structure interaction solution would be more precise for the solution process. 

Even with several questionable physical models for particular aspects of the 
whipping model, comparison with test data shows that the model provides reasonable 
results, particularly in early-time response. It is believed that the general approach of 
solving the equations of motion directly in the time domain was reasonable and necessary. 
It is in the modeling of each individual component of the forces on the whipping ship 
where future work should be concentrated. 
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% Subroutine IN.M 
% 

% This matiab subroutine is written as an INPUT FILE for all other MATLAB 
% subroutines. Data is loaded from formatted batch-type input files using 
% the MATLAB ’’load” command, and arrays and constants are generated for 
% the other MATLAB subroutines. 

% 

% MISCELLANEOUS/DEFAULT VALUES 

% 

% Number of hull stations or finite beam elements 
NE=19; 

% Number of hull sections or finite beam nodes 
NNODES=20: 

% Number of degrees of freedom (4 DOF at each node) 

NDOF=4*NNODES; 

% 

% 

% PROPERTIES OF EACH FINITE BEAM ELEMENT 
load element.dat 

% 

% Element number 
e=element(L:): 

% Cross-sectional area of element (in2) 

Ae=element(2.:); 

% Length of element (ft) 

Le=element(3,:); 

% Area effective in y-directed shear (in2) 

Asy=element(4,:); 

% Area effective in z-directed shear (in2) 

Asz=element(5,:); 

% Second moment of area about y-axis (in2ft2) 

Iav~ element(6.:); 

% Second moment of area about z-a.\is (in2ft2) 

Iaz=element(7.:): 

% Y-distance from the neutral axis to the outer fibet of the crown (ft) 
Yc=element(8.:): 

% Y-distance from the neutral axis to the outer fiber of the keel (ft) 
Yk=element(9.:); 

% Z-distance from the centerline to the outer fiber of the port hull (ft) 
Zp=element(10.:); 

% Z-distance from the centerline to the outer fiber of the stbd hull (ft) 
Zs=element(l 1.:); 

% 

% 

% PROPERTIES OF EACH HULL SECTION (NODE) 
load node.dat 
% 

% Node number 
nnode=node(l.:): 

% Length of each hull section (ft) 

Lh=node(2,:); 

% Weight of hull section (i.e. in air) (lton) 
vvh=node(3.:); 

% mass of hull section (lbf-sec2/ft=slug) 
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mh=wh*2240./32. 174; 

% Buoyancy of hull section (displacement of salt water) (lton) 
bh=node(4.:); 

% Weight moment of inertia of hull section about y-axis (lton-ft2) 
Imyh=node(5.:): 

% Mass moment of inenia of hull section about y-axis (slug-ft2) 
Imy t =Imyh*2240./32. 1 74; 

% Weight moment of inertia of hull section about z-axis (lton-ft2) 
Imzh=node(6,:); 

% Mass moment of inertia of hull section about z-axis (slug-ft2) 
Imz=Imzh*2240./32. 174; 

% Maximum section beam of main hull perp. to motion in z direction (ft) 
Bz=node(7.:); 

% Maximum section beam of main hull perp to motion in y direction (ft) 
B\-node(8.:); 

% Cross-sectional area of Appended hull (appended portion of each section) 
% perp. to motion in y direction (ft2) 

Aapv=node(9.:); 

% Cross-sectional area of Appended hull (appended portion of each section) 
% perp. to motion in z direction (ft2) 

Aapz=node(10.:); 

% Width of Appended hull perp. to motion in y direction (ft) 

Dapy=node(l 1.;); 

% Width of Appended hull perp. to motion in z direction (ft) 
Dapz=node(12,:); 

% Length of Appended hull perp. to motion in y direction (ft) 
Lapy=node(13,:); 

% Length of Appended hull perp. to motion in z direction (ft) 
Lapz=node(14.:): 

% 

% Hydrodynamic coefficients (non-dimensional) 

% 

% Added mass, inertia, and drag coefficients for each main hull section 
% for motion in vertical (y) and horizontal (z) directions 
Caly=node(15.:); 

Calz=node(16,:); 

Cmly=node(17,:); 

Cmlz=node(18.:); 

Cdl\ r =node(19 ;> :): 

Cdlz=node(20,:); 

% Added mass, inertia, and drag coefficients for appended hull sections 
% (appended portion of each hull section) for motion in y and z directions 
Ca2y=node(21.:); 

Ca2z=node(22.:); 

Cm2v=node(23.:); 

Cm2z=node(24.:); 

Cd2y-node(25.:). 

Cd2z=node(26.:); 

% 

% 

% MISCELLANEOUS INPUT 
load misc.dat 

% 

% Acceleration due to gravity (ft/sec2) 
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g=misc( 1); 

% Density of salt water (Ibf-sec2/ft4=slugs/ft3) 
rho=misc(2); 

% Young's modulus for steel (Ibf7in2) 

E=misc(3); 

% Poisson ratio for steel (non-dimensional) 

Nu=misc(4): 

% Yield stress for steel (Ibf7in2) 

Sig}^misc(5); 

% Shear modulus for steel (lbf/in2) 

G=misc(6); 

% Assumed material (hysteresis) modal viscous damping coefficients for 
% all low frequency (fundamental) flexural "dry " modes (non-dimensional) 
damp=misc(7); 

% Whipping mode to scale as initial condition for use in FR£E\VHIP.M 
MODE=misc(8); 

% Switch to determine whether or not "wet" whipping is used in FREE WHIP.M 
WET=misc(9); 

% Node numbers to plot displacements, velocities.accelerations 
NODEl=misc(10); 

NODE2=misc(ll); 

NODE3=misc(12); 

NODE4=misc(13); 

% Element numbers to plot strains 
ELEMENTl=misc(14); 

ELEMENT2=misc( 15); 

ELEMENT3=misc( 16); 

% 

% Input parameters for Newmark time integration algorithms 

% 

% time step (sec) 
dt=misc(17); 

% gamma parameter (non-dimensional) 
gammap=misc(18); 

% beta parameter (non-dimensional) 
betap=misc(19); 

% time limit (sec) 
tmax=misc(20); 

% vector norm tolerance 
tol=misc(21); 

% 

% Input for internal vibration absorbers 

% 

% "Reasonable" weight of an absorber or internal structure (lton) 

\va=misc(22): 

% mass of absorber or internal structure (slug) 
ma=wa*2240/32. 174; 

% Node number where absorber is "attached" to the hull 
nabs=misc(23): 

% 

% Input for internal sloshing tank 

% 

% "Reasonable" weight of liquid in a (large) internal tank (lton) 
wl=misc(24); 
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% mass of liquid (slug) 

Ml=\vl*2240./32. 174; 

% "Reasonable” aspect ratio of tank depth/width (h/a) (non-dimensional) 
ha=misc(25); 

% Node number where tank is "attached" to the hull 
ntank=misc(26); 

% 

% Input for explosion bubble dynamics 

% 

% Initial depth of charge (ft) 

Dchg=misc(27); 

% Charge w eight of TNT (lb) 

Wchg=misc(28); 

% Initial location of charge relative to ship hull (in hull coordinate system) 

xchg=misc(29); 

ychg=misc(30); 

zchg=misc(31); 

% Initial depth of the ship (ft) 

Dship=misc(32); 
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% SUBROUTINE HULL.M 

% 

%This subroutine performs the following 



Calculates the global mass and stiffness matrices (K and M) 

for the 3-D beam finite element formulation (ignoring axial 
and torsional effects). 

Solves the generalized eigenvalue problem to obtain the 
"dry mode" eigenvalues (natural freqs squared) and 
corresponding eigenvectors ("dry mode" shapes). 

Calculates the Rayleigh coefficients for critical hull 

(structural and hull hysteresis) damping defined by 
the 1st and 3rd lowest fundamental flexural modes 

Calculates the equivalent viscous hull damping matrix for the 
calculated Rayleigh coefficients. 

Calculates vectors for modal mass, modal damping, and 

modal stiffness (for N modes) corresponding to the "dry 
mode" shape vectors. These are used for optimization 
of internal absorber and liquid sloshing dampers (tanks) 
in other subroutines 

Prov ides for validation against the "exact" solution given by 
Rao (chapter 8) for a continuous uniform free-free beam. 
The natural frequencies of the finite element solution 
can be compared to the natural frequencies of the 
"exact" solution. This may be disabled after validation 
by placing % before each calculation. 



% 1 

% 

% 

% 2 . 

% 

% 

% 3. 

% 

% 

% 4 . 

% 

% 5. 

% 

% 

% 

% 

% 6 . 

% 

% 

% 

% 

% 

% 

%(D 

% Formulate the global (hull) stiffness and mass matrices (K and M). 

% 

% Initialize K and M 
K=zeros(NDOF); 

M=zeros(NDOF): 

% 

% Determine the global stiffness matrix (Ibf/ft for translation or or lbf-ft for rotation) 

% 

% Determine the element stiffness matrices (KE) 

% 

% Step through each of the 19 elastic beam elements 
for e=l:NE 



% 

Pv=(I2*E*Iaz(e))/(G*Asy(e)*Le(e) A 2); 

Pz=(12*E*Iay(e))/(G*Asz(e)*Le(e) A 2); 
bz=E*Iaz(e)/(( l+Pv)*Le(e) A 3); 
by-E*Iay(e)/(( I+Pz)*Le(e) A 3): 

K1 l=[12*bz 0 0 6*bz*Le(e); 0 12*by -6*by*Le(e) 0; 

0 -6*by*Le(e) (4+Pz)*by*Le(e) A 2 0; 6*bz*Le(e) 0 0 (4+Pv)*bz*Le(e) A 2|; 
K2I=[-12*bz 0 0 -6*bz*Le(e); 0 -12*by 6*by*Le(e) 0: 

0 -6*by*Le(e) (2-Pz)*by*Le(e) A 2 0; 6*bz*Le(e) 0 0 (2-Py)*bz*Le(e) A 2]; 
K12=K2F; 

K22=(I2*bz 0 0 -6*bz*Le(e): 0 12*by 6*by*Le(e) 0: 

0 6*by*Le(e) (4+Pz)*by*Le(e) A 2 0; -6*bz*Le(e) 0 0 (4+Py)*bz*Le(e) A 2]; 



KE=[K11 K12; K21 K221; 

% 



113 



% Transform element stiffness matrices into their 
% global coordinate matrices 
for ie=l:8 

for je=l:8 

i=ie+4*(e-l): 

j=je+4*(e-l); 

K(ij)=K(i.j)+KE(ieje): 

end 

end 

end 

% 

% Once each of the elastic beam elements has been considered and added. 

% the global (hull) stiffness matrix is K. 

% 

% Determine the global mass matrix (lbf-sec2/ft=slug for translation or slug-ft2 for rotation) 

% 

% Step through each of the 20 nodes 

% 

for n=l:NNODES 

% Step through each of the 4 dof at each node and define the NODE mass matrix 
MN(l,l)=mh(n); 

MN(2.2)=mh(n); 

MN(3,3)=Imy(n); 

MN(4.4)=Imz(n); 

% Transform node mass matrices for each node into their global coordinate matrices 
for in=l:4 

for jn=l:4 

i=in+4*(n-l); 

j=jn+4*(n-l): 

M(ij)=M(ij)+MN(in.jn); 

end 

end 

end 

% 

% Once each of the nodes has been considered, the global 
% mass matrix is M 

% 

% ( 2 ) 

% Solve the eigenvalue problem to obtain the "dry mode" natural frequencies (square roots of 
% eigenv alues) and corresponding "dry mode" shape vectors (eigenvectors). 

% 

% *** It should be noted that because the hull has 4 rigid body modes, 

% the lowest 4 eigenvalues will be zero (zero natural frequency). 

% These correspond to the lowest 4 eigenvectors which are the 
% RIGID BODY mode shapes. Thus, the lowest FLEXURAL mode 
% w ould be mode n=5. 

% 

% "Phil" is a (NxN) matrix whose columns are "dry mode" shape vectors. 

% "Wl" is a diagonal (NxN) matrix whose diagonal elements are (square of the natural frequencies). 
% (Note: W is NOT necessarilv defined in ascending order). 

% 

|Phil.Wl]=eig(K.M); 

% 

% Define a VECTOR of natural frequencies: 
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for r=l NDOF 

w l(r)=sqrt(Wl(r.r)); 
end 

% 

% Sort natural frequencies and define a vector of natural frequencies giv en in ascending order: 
|w2.i|=sort(wl); 

% where "i" is a vector of indices such that vv2=vv l(i). 

% 

% Sort the matrix of eigenvectors (mode shapes) as columns in an order corresponding to ascending 
% natural frequencies (defined by index vector "i"): 
for r=l:NDOF 
s=i(r); 

Phi2(:,r)=Phil(:,s); 

end 



% 



% Plot the 2 rigid body mode shapes and the first 4 flexural mode shapes corresponding to motion in the 
% horizontal plane (for example): 

% Step through the first 12 modes: 

%for m= 1:12 
% phi=Phi2(:.m). 

% for n=l:NNODES 

% nd(n)=n: 

% ny=l+4*(n-l): 

% nz=2+4*(n-l); 

% phiy(n)=phi(ny); 

% phiz(n)=phi(nz); 

% end 

% % For rigid body modes (m<=4): 

% if m<=4 



% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 



figure) 1) 
if m<=2, rbm=l; 
elseif m<=4. rbm=2: 
end 

if phiz( 1 ) — =0 

subplot(2,l.rbm), plot(nd.phiz) 
axis([l NNODES -1.0 1.0)) 

title('Dry mode shape of rigid body mode in horizontal (v=0) plane') 
end 

% For flexural modes (m>4): 
else 

figure(2) 
mcount=(m-4); 
if mcount<=2 

flexmode=l; 
elseif mcount<=4 

flexmode=2: 
elseif mcount<=6 

flexmode=3: 
else. flexmode=4; 
end 

if phiz( 1 )~=0 

subplot(4. 1 .flexmode), plot(nd.phiz) 
axis([l NNODES -1.0 1.0)) 

title(['Dry mode shape of flexural mode '.int2str( flexmode).' in ... 
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end 



horizontal (y=0) plane']) 



% 

% 

% 



end 



% end 

% Print "dry" mode shapes 
% figure(l) 

% set(gcf,'PaperPosition',[l.(). 4.5, 6.5, 4.0|) 

% print 

% figure(2) 

% set(gcf,’PaperPosition'.[l.(). 2.0, 6.5, 8.5]) 

% print 

% 

% 

% (3) 

% Calculate the Rayleigh coefficients defined by the natural frequencies of the 1st and 3rd flexural modes 
% (1st and 2nd flexural mode for one of the planes of motion - horizontal or vertical). Ignoring the 4 
% rigid bodv modes, these would be correspond to modes n=5 and n=7. 

% 

% Define the "dry mode" natural frequencies for modes 5 and 7: 

\v5=w2(5); 

vv7=w2(7); 

% 



% Find the Rayleigh coefficients using equation (3-16): 

w57=[vv7 -w5; -l/\\7 l/\v5]; 

damping=[damp; damp], 

pre=2*(vv5*vv7)/(\v7 A 2-\v5 A 2); 

coefF=pre*\v57*damping; 

alpha=coeff(l); 

beta=coeff(2); 

% 



% 



% (4) 

% Calculate the equivalent viscous damping matrix for hull damping using the calculated Rayleigh 
% coefficients: 

% 



C=alpha*M+beta*K; 

% 

% 

%(5) 

% Define the "dry " mode shape vectors and natural frequencies and calculate vectors for modal mass 
% (Mn), modal damping (Cn), and modal stiffness (Kn) corresponding to the "dry mode" mode 
% shape vectors: 

% 



for n=l:NDOF 

phin=Phi2(:,n); 

phint=phin'; 

Phi(:,n)=phin; 

Mn(n)=phint*M*phin; 

Cn(n)=phint*C*phin; 

Kn(n)=phint*K*phin; 

w(n)=w2(n); 



% "dry " mode shape vectors 
% modal mass v ector 
% modal damping vector 
% modal stiffness vector 
% "dry " natural frequencies 



end 

% 

% (6) 



116 



% Validation of a sample finite element formulation (uniform beam) against the "exact" solution, for a 
% defined axisv mmetric (about x-axis) beam arrangement 
% 

% The "exact" solution is given by Rao (Chapter 8) for a uniform, continuous 
% free-free beam with axisymmetry about the x-axis. This "exact" solution 
% does not include shear effects, so the approximate solution should 
% provide natural frequencies which are slightly lower than the "exact". 

% 

% Natural frequencies for the "exact" solution for the first 3 transverse flexural modes (modes 1.2.3): 

% 

% Bnl 1=4.73; 

% Bnl2=7.85; 

% Bnl3=l 1.0; 

% rho=mh(l)/Lh(l); 

% Sqroot=sqrt((E*Iay(l))/(rho*(Lh( 1)*NE) A 4)); 

% omega 1 =(B nl 1 A 2 ) * Sqroot ; 

% omega2=(Bnl2 A 2)*Sqroot; 

% omega3=(Bnl3 A 2)*Sqroot; 

% 

% Compare "exact" solution to the (approximate) numerical solution. 

% For approximate solution, rigid body modes are ignored and the 1st. 2nd 
% and 3rd flexural modes for transverse displacement (in each plane (y=0.z=())) 

% are considered. It should be noted that because of axisymmetry, natural 
% frequencies and mode shapes in each plane should be the same 
% (i.e. repeated eigenvalues). 

% 

% wl=w(5); 

% w2=w(7); 

% \v3=w(9); 

% Compare 
% omega 1 

% wl 

% omega2 

% \v2 

% omega 3 

% vv3 



1 . Calculates the fluid added mass matrix for the submarine 

using a Morison equation approach (i.e. strip theory) This 
approach includes added mass effects in the translational 
degrees of freedom onh (rotational fluid added mass effects 
are not included) 

2. Adds the fluid added mass matrix to the hull mass matrix (subroutine 

HULL.M) and solves the generalized eigenvalue problem to obtain 
the "wet mode" eigenvalues (natural frequencies squared) and 
corresponding "wet mode" eigenvectors ("wet mode" shapes). 

3. Calculates vectors for "wet" modal mass, modal damping, and modal 

stiffness (for N modes) corresponding to the "wet mode" shape 
vectors. These are used for optimization of internal absorber 
and liquid sloshing dampers (tanks) in other subroutines. 



% SUBROUTINE WETMODE. M (called as a MATLAB function) 

% 

% This subroutine performs the following: 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

%(D 

% Calculate the fluid added mass matrix (lbf-sec2/ft=slug for translation). 

% 

% Initialize the fluid added mass matrix 
Ma=zeros(NDOF) ; 

% 

% Step through each of the 20 hull sections (nodes) and define the NODE fluid added mass matrices 
for n= UNNODES 

% Calculate the fluid added mass associated with the main hull (Mai) 

% and the appended hull (Ma2) 

% 

% Step through each of the 4 dof at each node and define the NODE 
% fluid added mass matrices for the main hull and appended hull 
MNa 1( 1 . 1 )=rho*(pi*By(n) A 2/4)*Ca 1 v(n)*Lh(n); 
MNa2(l,l)=rho*Aapy(n)*Ca2y(n)*Lapv(n); 
MNal(2.2)=rho*(pi 5,e Bz(n) A 2/4)*Calz(n)*Lh(n); 
MNa2(2,2)=rho*Aapz(n)*Ca2z(n)*Lapz(n); 

MNal(3,3)=0; 

MNa2(3,3)=0; 

MNal(4.4)=0; 

MNa2(4.4)=0; 

% Transform NODE fluid added matrices into their global coordinate matrices 
for in=l:4 

for jn=l:4 

i=in+4*(n-l); 

j=jn+4*(n-l); 

Ma(i.j)=Ma(ij)+MNal(in,jn)+MNa2(injn); 



end 



end 



end 

% 

% Once each of the nodes has been considered, the global fluid added mass matrix is Ma 

% 

%( 2 ) 

% Add the fluid added mass matrix to the hull mass matrix and solve the generalized eigenvalue 
% problem to obtain the "wet mode" eigenvalues (natural frequencies squared) and corresponding 



% "wet mode” eigenvectors ("wet mode" shapes). 

% 

% Total mass matrix (Mt): 

Mt=M+Ma; 

% 

% "Phil" is a (NxN) matrix whose columns are "wet mode" shape vectors. 

% "Wl" is a diagonal (NxN) matrix whose diagonal elements are (square of the "wet" natural 
% frequencies). (Note: W is NOT necessarily defined in ascending order). 

% 

|Phil.Wl]=eig(K.Mt): 

% 

% Define a VECTOR of natural frequencies: 
for r=l:NDOF 

vvl(r)=sqrt(Wl(r,r)); 

end 

% 

% Sort natural frequencies and define a vector of "wet" natural frequencies given in ascending order. 
[w2,i]=sort(wl); 

% where "i" is a vector of indices such that w2=wl(i). 

% 

% Sort the matrix of eigenv ectors ("wet" mode shapes) as columns in an order corresponding to ascending 
% "wet" natural frequencies (defined by index vector "i"): 
for r=l :NDOF 
s=i(r); 

Phi2(:.r)=Phil(:,s); 

end 

% 

% Plot the 2 rigid body "wet" mode shapes and the first 4 flexural "wet" mode shapes corresponding to 
% motion in the horizontal plane (for example): 

% Step through the first 12 modes: 

%for m=l:12 
% phi=Phi2(:.m); 

% for n=l:NNODES 

% nd(n)=n; 

% ny=l+4*(n-l); 

% nz=2+4*(n-l); 

% phiy(n)=phi(ny); 

% phiz(n)=phi(nz); 

% end 

% % For ngid body modes (m<=4): 

% 

% if m<=4 

% figure(l) 

% if m<=2, rbm=l; 

% elseif m<=4, rbm=2; 

% end 

% if any(phiz)~=0 

% subplot(2.1,rbm), plot(nd.phiz) 

% axis([l NNODES -1.0 1.0]) 

% title('Wet mode shape of rigid body mode in horizontal (y^) plane') 

% end 

% % For flexural modes (m>4): 

% else 

% figure(2) 
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% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 



mcount=(m-4); 
if mcount<=2 

flexmode=l, 
elseif mcount<=4 

flcxmode=2: 
elseif mcount<=6 

flexmode=3; 

else 

flexmode=4; 

end 

if phiz(l)— 0 

subplot(4,l, flexmode), plot(nd.phiz) 
axis([l NNODES -1,0 1,0|) 

title(['Wet mode shape of flexural mode \int2str(flexmode).’ in 
...horizontal (>-()) plane' |) 
end 
end 



% end 

% Print wet mode shapes 
% figure(l) 

% set(gcf,'PaperPosition'.[1.0. 4.5, 6.5. 4.0]) 

% print 

% figure(2) 

% set(gcf,'PaperPosition'.[1.0, 2.0, 6.5. 8.5]) 

% print 

% 



% ( 3 ) 

% Define the "wet" mode shape vectors and natural frequencies the first and calculate vectors for "wet" 
% modal mass (Mnw), modal damping (Cnvv), and modal stiffness (Kmv) corresponding to the "wet" 
% mode shape vectors using equations (2-13): 



% 



for n=l:NDOF 

phin=Phi2(:.n); 

phint=phin'; 

Phiw(:,n)=phin; 

Mnw(n)=phint*Mt*phin; 

Cnw(n)=phint*C*phin; 

Knw(n)=phint*K*phin; 

ww(n)=w2(n); 

end 

% 

% 

% 



% "wet" mode shape vectors 
% "wet" modal mass vector 
% "wet" modal damping vector 
% "wet" modal stiffness vector 
% "wet" natural frequencies 
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% FREE WHIP. M 

% 

% This subroutine performs the following 

% Calculates the response (nodal point displacements, velocities. 

% and accelerations) of the submarine dunng a 

% free vibration with initial displacements given by a 

% scaled mode shape vector for selected mode shapes. 

% Response is calculated using the New mark time integration 
% method, with iteration conducted at each time step using a 

% Newton-Raphson method. 

% Utilize the following: 

% - K.C.M.Phi calculated in HULL.M 

% - Hydrodynamic coefficients provided in 1N.M 

% - Parameters for Newmark algorithm provided in IN.M 

% 

% Define nodes to use for response calculations 
N(l)=NODEl; 

N(2)=NODE2; 

N(3)=NODE3. 

N(4)=NODE4; 

% Define elements to use for axial strain calculations 
E11=ELEMENT1; 

E12=ELEMENT2; 

E13=ELEMENT3; 

% 

% Set initial and final times for Newmark integration, and limit on number 

% of Newton-Raphson iterations at each time step 

t=0; 

tf=tmax; 

jlimit=50; 

% 

% Initialize the loading forces on the hull 
% Loading forces related to displacement 
Fdisp=zeros( 1 ,NDOF)’; 

% Loading forces related to velocity 
Fvel=zeros( LNDOF)'; 

% Loading forces related to acceleration 
Facc=zeros( LNDOF)’; 

% Calculate (constant) static forces (i.e. hydrostatic/buoyancy and weight) 
Fgra\ -zeros( 1 .NDOF)'; 

Fhs=zeros( 1 .NDOF)’; 
for n=l:NNODES 

ny=(n-l)*4+l; 

% Hydrostatic/buoyancy' force (Ibf) 

Fhs(nv)=bh(n)*2240; 

% Static force of gravity /w eight (lbf) 

Fgrav(ny)=-wh(n)*2240; 

end 

% Total loading forces 
Floadt=Fdisp+Fvel+Facc+Fhs+Fgrav; 

% 

% Define initial conditions (displacements/velocities/accelerations) 

% given by a scaled mode shape vector (for free w hipping) 
if MODE-=l. xt=(0.()2)*sum(Lh) 5<e Phi(:.6); 
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else. xt=(0.02)*sum(Lh)*Phi(:,8); 

end 

% 

xdt=zeros( LNDOF)’; 

Minv-inv(M); 

xddt=(Minv*(Floadt-C*xdt-K*xt)); 

% 

% Perform the Newmark time integration 

% 

% Step through time 

% 

% Define constants to be used in Neuton-Raphson iterations 

cl=(l/(betap*dt A 2)); 

c2=(l/(betap*dt)); 

c3=(l/(2*betap)); 

c4=(gammap/(betap*dt)); 

c5=gammap/betap; 

c6=( l-gammap/(2*betap)); 

c7=(l/(2*betap)-l); 

% 

dx=zeros( LNDOF)'; 

% 

for i=l:300; 

t=t+dt 

% Set initial displacments.velocities.accelerations, and loading forces for 

% the current time step (equal to the final for the last time step) 

xj=xt; 

xdj=xdt; 

xddj=xddt; 

Floadj=Floadt; 

% 

% Perform Newton-Raphson iterations at current time step 
for j= 1 :jlimit 

% 

% Calculate loading forces for the current iteration 
for n=l:NNODES 

% Loading forces related to velocity 
nv=4*(n-l)+l; 
nz=4*(n-l)+2; 

Fvel(ny)=-0.5*rho*Bv(n)*Cdly(n)*xdj(ny)*abs(xdj(ny))*... 
...Lh(n)-0.5*rho*Dapy(n)*Cd2y(n)*xdj(ny)*abs(xdj(ny))*Lapy(n); 
Fvel(nz)=-0.5*rho*Bz(n)*Cdlz(n)*xdj(nz)*abs(xdj(nz))*... 
...Lh(n)-0.5 3 * c rho*Dapz(n)*Cd2z(n) 5 * c xdj(nz)*abs(xdj(nz))*Lapz(n); 
% Loading forces related to acceleration 

Facc(ny)=-rho*(pi*By(n) A 2/4)*Caly(n)*Lh(n)*... 

...xddj(ny)-rho*Aapy(n)*Ca2y(n)*Lapy(n)*xddj(ny); 

Facc(nz)=-rho*(pi*Bz(n) A 2/4)*Calz(n)*Lh(n) 3 * c ... 

. f .xddj(nz)-rho*Aapz(n)*Ca2z(n)*Lapz(n)*xddj(nz); 

end 

% 

% Total loading force for the current iteration 
Floadj=Fdisp+Fvel+Facc+Fhs+Fgrav; 

% 

% Calculate Jacobian matrices for each of the loading forces 
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% (square diagonal matrices for this formulation) 

dFdxdj l=zeros( 1 .NDOF)'; 

dFdxddj l=zeros( 1 .NDOF)'; 

dFdxj=zeros(NDOF); 

dFdxdj=zeros(NDOF); 

dFdxddj=zeros(NDOF); 

dxdj=zeros( 1 .NDOF)'; 

dxddj=zeros( l,NDOF)’; 

for n=l:NNODES 

n\^4*(n-l)+l; 

nz=4*(n-l)+2; 

dxdj(ny)=0.01: 

dxdj(nz)=0.01; 

dxddj(ny)=0.01: 

dxddj(nz)=0.01; 

dFdxdj l(ny)= .. 

...(-0.5*rho*By(n)*Cdly(n)*(xdj(ny)+dxdj(ny))*abs(xdj(ny)+dxdj(nv))* 

...Lh(n)- 

0.5*rho*Dapy(n)*Cd2y(n)*(xdj(ny)+dxdj(ny))*abs(xdj(nv)+dxdj(ny))*... 
...Lapy(n)-Fvel(ny))/dxdj(ny); 
dFdxdj l(nz)= .. 

...(o*rho*Bz(n)*Cdlz(n)*(xdj(nz)+dxdj(nz))*abs(xdj(nz)+dxdj(nz))*.. 

...Lh(n)- 

0.5*rho*Dapz(n)*Cd2z(n)*(xdj(nz)+dxdj(nz))*abs(xdj(nz)+dxdj(nz))*... 

. . . Lapz(n)-Fvcl(nz))/dxdj(nz); 

dFdxddj l(ny)=(-rho*(pi*By(n) A 2/4)*Caly(n)*Lh(n)*... 
...(xddj(ny)+dxddj(ny))-rho*Aapy(n)*Ca2y(n)*Lapy(n)*... 
...(xddj(ny)+dxddj(ny))-Facc(ny))/dxddj(ny); 
dFdxddj l(nz)=(-rho*(pi*Bz(n) A 2/4)*Calz(n)*Lh(n)*... 
...(xddj(nz)+dxddj(nz))-rho*Aapz(n)*Ca2z(n)*Lapz(n)*... 
...(xddj(nz)+dxddj(nz))-Facc(nz))/dxddj(nz); 
end 

for n=l:NDOF 

dFdxdj(n.n)=dFdxdjl(n); 
dFdxddj(n,n)=dFdxddj l(n); 

end 

% 

% Calculate residuals 

flj=M*(cl%\j-xt)-c2^xdt-c3*\ddt)+C*(c4*(xj-xt)-c5*xdt+c6*xddt*dt)+... 

...K*(xj-xt)-Floadj+Floadt; 

f2j=xdj-c4*(xj-xt)+(c5-l)*xdt-c6*dt*xddt; 

f3j=xddj-c 1 *(xj-xt)+c2*xdt+c7*xddt; 

% 

% Solve for the displacement correction (dx) 

% 

A l=(c 1 *M+c4*C+K-c 1 * dFdxddj <4 * dF dxdj -dFdxj ) ; 

Alinv=inv(Al); 

A2=-fl j -dFdxdj * f2j-dFdxddj * Oj ; 
dx=(Alinv*A2); 

% 

% Calculate the next iteration's displacements, velocities, and accelerations 

% 

xj=xj+dx; 

xdj=xdj+(c4*dx)-f2j; 
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xddj=xddj+(cl*dx)-f3j; 

% 

% Check to see if all Norms of residuals are less than tolerance for 
% the current iteration 
% (If so. go to the next time step) 

Normflj=norm(flj); 

Normf2j=norm(f2j); 

Normf3j=norm(f3j); 

if Normfl j<=tol & Normf2j<=tol & Normf3j<=tol 

% 

% Save current time, displacements, velocities, 

% and accelerations (at specified nodes) for later plotting 
time(i)=t; 

% Displacements/velocities/accelerations 
for n=l:4 

Dh(n)=(N(n)- 1 )*4+2; 

end 

Xl(i)=xj(Dh( 1)); 

Xdl(i)=xdj(Dh(l)); 

Xdd 1 ( i)=xddj ( D h( 1 ) ) ; 

X2(i)=xj(Dh(2)); 

Xd2(i)=xdj(Dh(2)); 

Xdd2(i)=xddj(Dh(2)); 

X3(i)=xj(Dh(3)); 

Xd3(i)=xdj(Dh(3)); 

Xdd3(i)=xddj(Dh(3)); 

X4(i)=xj(Dh(4)); 

Xd4(i)=xdj(Dh(4)); 

Xdd4(i)=xddj(Dh(4)); 

% 

Dhl9=(19-l)*4+2; 

Dhl8=(18-l)*4+2; 

Fvel 1 8(i)=Fvel(Dh 18); 

Faccl8(i)=Facc(Dhl8); 

Floadl8(i)=Floadj(Dhl8); 

Fvel 1 9(i)=Fvel(Dh 1 9); 

Face 1 9(i)=Facc(Dh 1 9); 

Fload 1 9(i)=Floadj(Dh 1 9); 

% 

% Save the current time step displacements, velocities, and 

% accelerations, and forces for use in the next time step 

xt=xj; 

xdt=xdj; 

xddt=xddj; 

Floadt=Floadj; 
break, end 
end 

if j==jlimit, break, end 
end 

if j==jlimit, break, end 
if t>=tf, break, end 
end 
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% Subroutine BUBBLE. M 

% 

% This subroutine performs the following: 

% 1. Computes the time histories of the explosion bubble radius (rad) 

% and bubble depth (depth) by integration of the differential 

% equations governing bubble dynamics. 

% Computation is earned out for specified bubble initial conditions 

% (i.e. initial depth of the charge (Dchg) and weight of the 

% charge (Wchg)). 

% 

% 

function[TMBR,BD,BRD.BDD.3RDD.BDDD.DEL.btime.brad,bdpth.bveLLl.Tl ) 
...bubble(Dchg.Wchg.tmax) 

% 

% Define constants 

% 

% Drag coeff. for bubble 
Cd=2.25; 

% Initial bubble head (ft) 

D0=Dchg+33; 

% Adiabatic gas exponent for TNT 
g 1 = 1 -25; 

% Constant for TNT 
kl=(0.0552)*D() A 0.25; 

% Length scale constant 
L 1 = 13. 67*(Wchg/D0) A ( 1/3); 

% Time scale constant 
Tl=2.94*(Wchg A (l/3)/D0 A (5/6)); 

% Misc. constants 
cl=33/Ll; 
c2=D0/Ll; 
c3=(gl-l)*kl; 
c4=3*gl + l; 

% 

% Define initial and final conditions 

% 

% Non-dim. bubble radius 
xO=k 1 A (4/3)*( 1 +(4*k 1 A 4)/3); 

% Non-dim. bubble depth 
y0=c2; 

% Rate of change of bubble radius 
xd0=0; 

% Rate of change of bubble depth 
yd0=0; 

% Non-dim. start time (i.e. 0 seconds) 

T0=0.(); 

% Non-dim. end time (i.e. tmax seconds) 

Tf=tmax/Tl; 

% Non-dimensional time step 
dT=( 5.0* xO A 4 )/T f; 

% Integration time step parameter 
tsp=10.; 

% 

% Perform Runge-Kutta time iteration 



125 



T=TO; 

x=xO; 

v=yO; 

xd=xdO; 

yd=ydO; 

% 

M=3000; % defined maximum number of time steps (for control) 

% 

for m=l:M 

% Call the differential eqn. for bubble radius 
[x,xd]=radfn(dT,x,xd,y.yd.c 1 ,c2.c3.c4.Cd); 

xdd=-1.5/(l-x/(2*(y-cl)))*((xd A 2/x)*(I-2^\/(3*(y-cl)))-yd A 2/(6*x)+y/(x*c2)-c3/(x A c4)+... 

..Xx/(4*(y-cl) A 2))*(xd*yd/3-x/c2+Cd*yd A 2/4)); 

% Call the differential eqn. for bubble depth 
(\ r ,yd]=dpthfn(dT.x,xd,y,yd,xdd.cl.c2.c3.c4.Cd): 

ydd=-3*(l/c2+xd*yd/x-(Cd*yd A 2)/(4^\)+(x/(4*(y-cl) A 2))*(3*xd A 2+x*xdd)); 

% 

% Save the current bubble characteristics for plotting and later use 

% 

TM(m)=T; % non-dim. time 

BR(m)=x; % non-dim. bubble radius 

BD(m)=v; % non-dim. bubble depth 

BRD(m)=xd; % non-dim. rate of change of bubble radius 

BDD(m)=yd; % non-dim. rate of change of bubble depth (bubble velocity) 

BRDD(m)=xdd; 

BDDD(m)=ydd; 

DEL(m)=(y-cl); % non-dim. suface pressure 
% 

btime(m)=T*Tl; % dimensional time 
brad(m)=x*Ll; % dimensional bubble radius 
bdpth(m)=(y-cl)*Ll; % dimensional bubble depth 
bvel(m)=vd*Ll/Tl; % dimensional bubble velocity 

% 

% Go to next time increment 
x3=x A 3; 

x3d=3.0*(x A 2)*xd; 

x3dd=3.0*x*(x*xdd+2.0*xd A 2); 

dT=2.0/(abs(x3dd)*3.()*tsp); 

if dT<0.001, dT=0.0()l; 
elseif dT>0.025. dT=0.025; 
end 

T=T+dT: 

if T>Tf. break 
end 
end 

% 

% Plot the results 

% 

figured ) 

subpIot(3. 1. 1). plot(btime.brad) 
title('Bubble radius vs. time') 
xlabel('Time (sec)') 
ylabel('Bubble radius (ft)') 
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subplot(3. 1.2). plot(btime.-bdpth) 
title('Bubble depth vs. time') 

\label('Time (sec)') 
ylabel('Bubble depth (- fit)') 
subplot(3.1.3). plotfbtime.-bvel) 
title('Upuard bubble velocity vs. time') 
xlabelf'Time (sec)’) 
ylabcU'Upward bubble velocity ( ft/s)') 
set(gcf.'PaperPosition'.| 1.0. 2.0. 6.5, 8.51) 
% print 
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% RADFN.M 

% 

% Function M-file for calculation of the bubble radius 
% (i.e. contains the bubble radius governing equation for use in BUBBLE. M) 

% 

function [x.xd]=radfn(dT.x.xd.y.yd.c 1 .c2.c3.c4.Cd) 

akl=dT*rfunc(x.xd,y.yd,cl.c2.c3.c4,Cd). 

ak2=dT*rfunc(x+dT*xd/2.0.xd+akl/2.0.y.yd,cl.c2.c3.c4.Cd): 

ak3=dT*rfunc(x+dT*xd/2.0+dT*akl/4,0.xd+ak2/2.0.y.yd.cl.c2.c3.c4.Cd); 

ak4=dT*rfunc(x+dT*xd+dT*ak2/2.0.xd-t-ak3.y.yd.cl.c2.c3.c4.Cd). 

x=x+dT*xd+dT*(akl+ak2+ak3)/6.0; 

xd=xd+(akl+2.Q*ak2+2.0*ak3+ak4)/6.0; 



% RFUNC.M 
% 

function ak=rfunc(x.xd.y,yd,cl,c2.c3.c4.Cd) 

ak=-1.5/(l-x/(2*(y-cl)))*((xd A 2/x)*(l-2*x/(3*(y-cl)))-yd A 2/(6*x)+y/(x*c2)-c3/(x A c4)+(x/(4*(y- 

cl) A 2))*(xd*vd/3-x/c2+Cd*yd A 2/4)); 
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% DPTHFN.M 

% 

% Function M-file for calculation of the bubble depth 
% (i.e. contains the bubble depth governing equation for use in BUBBLE. M) 

% 

function [y.yd|=dpthfn(dT.x.xd.y.yd.xdd.cl.c2.c3.c4.Cd) 

akl=dT*dfunc(x.xd.y.yd.xdd.cl.c2.c3.c4.Cd); 

ak2=dT*dfunc(x.xd.y+dT*yd/2.0.yd+akl/2.0.xdd.cl.c2.c3.c4.Cd). 

ak3=dT*dfiinc(x.xd.y+dT*yd/2.0+dT*akl/4.0.yd+ak2/2.0.xdd.cl.c2.c3.c4.Cd); 

ak4=dT*dfunc(x,xd.y+dT*yd+dT*ak2/2.0,yd+ak3.xdd.cl.c2.c3.c4.Cd). 

\-=y+dT*yd+dT*(akl+ak2+ak3)/6.0; 

yd=> d+(ak 1 +2.0*ak2+2.0*ak3+ak4)/6.0; 



% DFUNC.M 
% 

% 



% 

function ak=dfunc(x.xd,y.yd.xdd.cl.c2.c3.c4.Cd) 

ak=-3*( l/c2+xd*yd/x-(Cd*yd A 2)/(4*x)+(x/(4*(y-cl) A 2))*(3*xd A 2+x*xdd)): 
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% BUB WHIP. M 

% 

% 

% This subroutine performs the following: 

% 1. Calculates the response (nodal point displacements, velocities, accelerations) 

% of the submarine during a vibration with initial displacements of zero 

% and subjected to fluid loading produced by a pulsating explosion bubble. 

% Response is calculated using the Newmark time integration algorithm, 

% with iteration conducted at each time step using a Nevvton-Raphson method. 

% 2. Calculates the local axial strains at specified hull locations using simple 

% beam theory (used for comparison with model test data). 

% 

% 

% Utilize the following: 

% -K.C.M calculated in HULL.M 
% -Hydrodynamic coefficients prov ided in 1N.M 
% -Parameters for Newmark algorithm provided in IN.M 
% -Bubble characteristics (radius, depth, etc.) calculated in BUBBLE. M 
% 

% Define nodes to use for response calculations 
N(l)=NODEl; 

N(2)=NODE2; 

N(3)=NODE3; 

N(4)=NODE4; 

% Define hull elements to use for axial strain calculations 
E11=ELEMENT1; 

E12=ELEMENT2; 

E13=ELEMENT3; 

% 

% Define the horizontal (hull) degree of freedom where internal structure/absorber acts 
nnabs=( nabs- 1 )*4+2 ; 

% 

% Define parameters for internal structure/absorber 

ka=ma*(vvw(8) A 2); 

mabs=zeros( LNDOF)’; 

kabs=zeros( 1 ,NDOF)'; 

mabs(nnabs)=ma; 

kabs(nnabs)=ka; 

% 

% Set the initial and final times (sec) for Newmark integration, and limit 

% on the number of Newton-Raphson iterations 

t=0; 

tf=tmax; 

jlimit=40; 

% 

% Initialize the loading forces on the hull 
% Loading forces related to displacement 
Fdisp=zeros( 1 ,NDOF)’; 

% Loading forces related to velocity 
Fvel=zeros( 1 ,NDOF)'; 

% Loading forces related to acceleration 
Facc=zeros( 1 .NDOF)’; 

% Calculate (constant) static forces (i.e. hydrostatic/buoyancy and weight) 

Fgrav=zeros( 1 ,NDOF)'; 
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Fhs=zeros( l.NDOF)'; 
for n=l:NNODES 

ny=(n-l)*4+l; 

% Hydrostatic/buoyancy force (lbf) 

Fhs(nv)=bh(n)*2240: 

% Static force of gravity/weight (lbO 
Fgrav(ny)=-vvh(n)*2240; 
end 

% Total initial loading forces 
Floadt=Fdisp+Fvel+Facc+Fhs+Fgrav; 

% 

% Define initial conditions (displacements/velocities/accelerations) 

% Hull 

xt=zeros( l.NDOF)': 
xdt=zeros( l.NDOF)'; 
xddt=zeros( 1 .NDOF)'; 

% Internal absorber 

xat=0: 

xadt=0; 

xaddt=0; 

% 

% Perform the Newmark time integration 

% 

% Step through time 

% 

% Define constants to be used in Nevvton-Raphson iterations 

cl=(l/(betap*dt A 2)): 

c2=(l/(betap*dt)); 

c3=(l/(2*betap)); 

c4=(gammap/(betap*dt)); 

c5=gammap/betap; 

c6=( l-gammap/(2*bctap)): 

c7=(l/(2*betap)-l). 

% 

dx=zeros( 1 .NDOF)'; 

% 

for i=l: 199; 

t=t+dt 

% 

% Set initial displacments.\elocities.accelerations. and loading forces for 

% the current time step (equal to the final for the last time step) 

xj=xt; 

xdj=xdt; 

xddj=xddt; 

Floadj=Floadt; 

xaj=xat; 

xadj=xadt; 

xadds=xaddt; 

% 

% Perform Nevvton-Raphson iterations at current time step 
for j=l:jlimit 

% 

xaddj=xadds: 

% Calculate current iteration absorber displacement/velocity 
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xaj=xat+xadt*dt+xaddt*(dt A 2/2)+betap*(xaddj-xaddt)*dt A 2:, 

xadj=xadt+xaddt*dt+gammap*(xaddj-xaddt)*dt; 

% 

% Calculate the loading forces for the current iteration 

% 

% Call function '’fluid” to provide free field fluid velocities and 
% accelerations at each node (at the current iteration node location) 

|uy,uz.uyd,uzd]=fluid(xj,xchg,vchg,zchg.Dship,uTM,BR,BD.BRD.BDD.BRDD.BDDD,. 

DEL.L 1 .T 1 .btime.bdpth.NDOF.NNODES.Lh); 

% 

for n=l:NNODES 

n\-4*(n-l)+l; 

nz=4*(n-l)+2; 

% Loading forces related to displacement 
Fdisp(nz)=kabs(nz)*(xaj-xj(nz)); 

% Loading forces related to velocity 

Fvel(ny)=0.5*rho*By(n)*Cdly(n)*(uy(n)-xdj(ny ))*... 

...abs(uy(n)-xdj(ny))*Lh(n)+()J*rho*Dapy(n)*Cd2y(n)*(uy(n)- 

...xdj(ny))*abs(uy(n)-xdj(ny))*Lapy(n); 

Fvel(nz)=0.5*rho*Bz(n)*Cdlz(n)*(uz(n)-xdj(nz))*... 

...abs(uz(n)-xdj(nz))*Lh(n)+0.5*rho*Dapz(n)*Cd2z(n)*(uz(n)- 

...xdj(nz))*abs(uz(n)-xdj(nz))*Lapz(n); 

% Loading forces related to acceleration 

Facc(ny)=rho*(pi + By(n) A 2/4)*Cmly(n)*Lh(n)*(uyd(n)- 

...xddj(ny))+rho*Aapv(n)*Cm2y(n)*Lapy(n)*(uyd(n)- 

...xddj(ny))+rho*(pi*By(n) A 2/4)*Lh(n)*xddj(ny)+rho*Aapy(n)*Lapy(n)*xddj(ny); 

Facc(nz)=rho*(pi*Bz(n) A 2/4)*Cmlz(n)*Lh(n)*(uzd(n)- 

...xddj(nz))+rho*Aapz(n)*Cm2z(n)*Lapz(n)*(uzd(n)- 

..xddj(nz))+rho*(pi*Bz(n) A 2/4)*Lh(n)*xddj(nz)+rho*Aapz(n)*Lapz(n)*xddj(nz); 

end 

% 

% Total loading force for the current iteration 
Floadj=Fdisp+Fvel+Facc+Fhs+Fgrav; 

% 

% Calculate Jacobian matrices for each of the loading forces 
% (square diagonal matrices for this formulation) 

% 

dFdxj 1 =zeros( 1 ,NDOF)'; 
dFdxdj 1 =zeros( 1 ,NDOF)’: 
dFdxddj 1 =zeros( 1 ,NDOF)'; 
dFdxj =zeros(NDOF) ; 
dFdxdj=zeros(NDOF); 
dFdxddj=zeros(NDOF); 
dxj=zeros( 1 ,NDOF)’; 

<±\dj=zeros( 1 ,NDOF)’; 
dxddj=zeros( LNDOF)'; 
for n=l :NNODES 

ny=4*(n-l)+l; 

nz=4*(n-l)+2; 

dxj(nz)=0.01; 

dxdj(ny)=0.01; 

dxdj(nz)=0.01; 

dxddj(ny)=0.01; 
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dxddj(nz)=0.01. 

dFdxj l(nz)=(kabs(nz)*(\aj-(xj(nz)+dxj(nz)))-Fdisp(nz))/d\j(nz); 
dFdxdj l(ny)=(0.5*rho*By(n)*Cdly(n)*(uy(nMxdj(ny)+dxdj(ny)))*abs(uy(n) 
...(\dj(ny)+d\dj(ny)))*Lh(n)H).5*rho*Dapy(n)*Cd2y(n)*(uy(n)- 
...(\dj(ny)+dxdj(ny)))*abs(uy(n)-(xdj(ny)+dxdj(ny)))*Lapy(n)- 
. Fvel(ny))/dxdj(ny); 

dFdxdj l(nz)=(0.5*rho*Bz(n)*Cdlz(n)*(uz(n)-(xaj(nz)+dxdj(nz)))*abs(uz(n) 
..Xxdj(nz)+dxdj(nz)))*Lh(n)+0.5*rho*Dapz(n)*Cd2z(n)*(uz(n)- 
...(xdj(nz)+d\dj(nz)))*abs(uz(n)-(xdj(nz)+dxdj(nz)))*Lapz(n)- 
...Fvel(nz))/dxdj(nz); 

dFd\ddjl(nv)=(rho*(pi*By(n) A 2/4)*Cmly(n)*Lh(n) 3 * : (uyd(n)- 

...(xddj(ny)+dxddj(ny)))+rho*Aapy(n)*Cm2y(n)*Lapy(n)*(uyd(n)- 

(xddj(ny)+dxddj(ny)))+rho*(pi*By(n) A 2/4)*Lh(n)*(xddj(ny)+dxddj(ny))+rho*Aapv(n)* 

...Lapy(n)*(xddj(ny)+dxddj(ny))-Facc(ny))/d\ddj(ny); 

dFdxddjl(nz)=(rho*(pi*Bz(n) A 2/4)*Cmlz4n)*Lh(n)*(uzd(n)- 

...(xddj(nz)+dxddj(nz)))+rho*Aapz(n)*Cm2z(n)*Lapz(n)*(uzd(n)- 

(xddj(nz4+dxddj(nz)))+rho*(pi*Bz4n) A 2/4)*Lh(n)*(xddj(nz)+dxddj(nz))+rho*Aapz(n)* 

...Lapz(n)*(xddj(nz)+dxddj(nz))-Facc(nz))/dxddj(nz); 

end 

for n=l:NDOF 

dFdxj(n,n)=dFdxjl(n); 
dFdxdj ( n , n )=dFdxdj 1 ( n) ; 
dFdxddj(n,n)=dFdxddjl(n); 
end 
% 

% Calculate residuals 

flj=M*(cl*(xj-xt)-c2*xdt-c3*xddt)+C*(c4*(xj-xt)-c5*xdt+c6*xddt*dt)+... 

. . . K*(xj-xt)-Floadj+Floadt; 

12j = xdj-c4*(xj-xt)+(c5-l)*xdt-c6*dt*xddt; 

f3j=xddj-cl%\j-xt)+c2*xdt+c7*xddt; 

% 

% Solve for the displacement correction (dx) 

% 

A l=(cl *M+c4*C+K-c 1 *dFdxddj-c4* dFdxdj -dFdxj); 

Alinv-inv(Al); 

A2=-f 1 j -dFdxdj * f2j -dFdxddj * Oj ; 
dx=(Alinv*A2); 

% 

% Calculate the next iteration's displacements, velocities, and accelerations 

% 

\j=\j+d\; 

xdj=xdj+(c4*dx)-f2j; 
xddj=xddj+(cl *dx)-f3j; 

% 

% Required acceleration of absorber 
if any(mabs)~=0, xadds=-Fdisp(nnabs)/mabs(nnabs); 
else, xadds=0; 
end 

% 

% Check to see if all Norms of residuals are less than tolerance for 
% the current iteration 
% (If so. go to the next time step) 

Normflj=norm(flj); 

Normf2j=norm(f2j); 
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Normf3j=norm(f3j); 

Normabs=norm( xadds-xaddj ) ; 

if Normflj<=tol & Normf2j<=tol & Normf3j<=tol & Normabs<=tol 

% 

% Calculate axial stress at specified hull stations 
% (using simple beam theory) 

% 

qy=zeros(l.NEy; 

qz=zeros(l.NE)'; 

V\-zeros( 1,NE)'; 

Vz=zeros(l,NEy; 

Myl=zeros( 1,NE)'; 

Mzl=zeros(l,NEy; 

My=zeros(l,NE)'; 

Mz=zeros(l,NE)’; 
for e=l:NE 

% Net vertical or horizontal force at element e 

ey=(e-l)*4+l; 

ez=(e-l)*4+2; 

qy(e)=Floadj(ey)-M(ey,ey)*xddj(ey); 

qz(e)=Floadj(ez)-M(ez,ez)*xddj(ez); 

% 

% Vertical or horizontal shear force at element e 
Vy(e)=sum(qy); 

Vz(e)=sum(qz); 

% 

% Vertical or horizontal bending moment at element e 
Myl(e)=Vy(e)*Le(e); 

Mzl(e)=Vz(e)*Le(e); 

My(e)=sum(Myl); % BM about z-axis 
Mz(e)=sum(Mzl); % BM about y-axis 

% 

% Axial stress at crown of element e (positive is tension) 
sigc(e)=-My(e)*Yc(e)/Iaz(e); 

% Axial stress at keel of element e 
sigk(e)=My(e)*Yk(e)/Iaz(e); 

% Axial stress at port-side fiber (port=pos z) (positive is tension) 
sigp(e)=-Mz(e) 5|c Zp(e)/Iay(e); 

% Axial stress at stbd-side fiber 
sigs(e)=Mz(e) *Zs(e)/Iay(e) ; 

% 

% Axial strains 
EPSc(e)=sigc(e)/E; 

EPSk(e)=sigk(e)/E; 

EPSp(e)=sigp(e)/E; 

EPSs(e)=sigs(e)/E; 

% 

end 

% Save current time, displacements and velocities 
% (at specified nodes) and axial strains 
% (at specified elements) for later plotting... also save 
% current horizontal fluid velocity and acceleration, and 
% (at nodes 10 and 20) and all calculated horizontal forces 
% (at nodes 10 and 20) for later plotting. 



134 



time(i)=t; 

% Vertical/horizontal dispIacementsA elocities 
% at specified nodes 
for n=l:4 

Dv(n)=(N(n)-l)*4+l. 

Dh(n)=(N(n)-l)*4+2: 

end 

xvl(i)=xj(Dv( 1 )): 

xhl(i)=xj(Dh( 1)): 

xd\i(i)=xdj(Dv(l)); 

xdhl(i)=xdj(Dh(l)): 

xv2(i)=xj(Dv(2)); 

xh2(i)=xj(Dh(2)): 

xdv2(i)=xdj(Dv(2)); 

xdh2(i)=xdj(Dh(2)); 

xv3(i)=xj(Dv(3 )); 

xh3(i)=xj(Dh(3)): 

xdv3 ( i)=xdj ( Dv( 3 )); 

xdh3(i)=xdj(Dh( 3 )): 

xv4(i)=xj(Dv(4)): 

xh4(i)=xj(Dh(4)): 

xdv4(i)=xdj(D\(4)); 

xdh4(i)=xdj(Dh(4)); 

% 

xabs(i)=xaj: 

xdabs(i)=xadj; 

xddabs(i)=xaddj; 

% 

% Horizontal hull/fluid motions and 

% external forces 

xddh2(i)=xddj(Dh(2)); 

xddh3 ( i)=xddj(D h( 3 )) ; 

uh3(i)=uz(N(3)); 

udh3(i)=uzd(N(3)): 

Fvelh3(i)=Fvel(Dh(3)); 

Facch3(i)=Facc(Dh(3)): 

Floadh3(i)=Floadj(Dh(3)): 

% 

xddh4(i)=xddj(Dh(4)); 

uh4(i)=uz(N(4)); 

udh4(i)=uzd(N(4)); 

Fvelh4(i )=Fvel(Dh(4»; 

Facch4(i)=Facc(Dh(4)); 

Floadh4(i )=Floadj ( D h( 4 ) ) ; 

% 

% Axial strains (keel, crown, port, stbd) at 

% at specified elements 

straink l(i)=EPSk(El 1 ): 

strainc 1 (i)=EPSc(El 1 ); 

strainpl(i)=EPSp(El 1 ): 

strains 1 (i)=EPSs(El 1 ); 

straink2(i)=EPSk(E12); 

strainc2(i)=EPSc(E12); 

strainp2(i)=EPSp(E12): 
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strains2(i)=EPSs(E12); 

straink3(i)=EPSk(E13); 

strainc3(i)=EPSc(E13); 

strainp3(i)=EPSp(E13); 

strains3(i)=EPSs(E13); 

% 

% Save the current time step displacements, velocities, and 

% accelerations, and forces for use in the next time step 

xt=xj. 

xdt=xdj: 

xddt=xddj; 

Floadt=Floadj; 
xat=xaj; 
xadt=xadj; 
xaddt=xaddj: 
break, end 
end 

ifj==jlimit. break, end 
end 

ifj==jli m it. break, end 
if t>=tf. break, end 
end 
% 
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% FLUID. M 

% 

% Calculates the free field fluid velocities and accelerations (at each node) 

% 

functionfuy.uz,uyd.uzd]=fluidvel(xj.xchg.vchg.zchg.Dship,t.TM.BR.BD.BRD,BDD... 

...,BRDD.BDDD.DEL.LLTl.btime,bdpth.NT)OF.NNODESXh) 

% 

% Non-dimensional bubble parameters 
TI=t/Tl; 

a=interp 1 (TM.BRTI); 
s=interp 1 (TMBRD.TI); 
sd=interp 1 (TM.BRDDTI); 
l=interp 1 (TM,BDD,TI); 
ld=interp 1 (TM.BDDD.TI); 
d=interp 1 (TM,DEL.TI): 

% Source strengths 

el=(Ll A 3)*(a A 2)*s/Tl; 

e2=-(Ll A 4/(2*Tl))*(a A 3*l); 

eld=(Ll A 3/Tl A 2)*(a A 2*sd+2*a*s A 2); 

e2d=-(Ll A 4/(2*Tl A 2))*(a A 3*ld+3*a A 2*l*s); 

% Vertical bubble velocity 
vm=-l*Ll/Tl; 

% Bubble depth 
ti=t; 

bdepth=interpl(btime.bdpth,ti); 

% 

for n=l:NNODES 

ny=(n-l)*4+l; 

nz=(n-l)*4+2; 

xnode=(n-l)*Lh(n); 

ynode=xj(ny); 

znode=xj(nz); 

% vr and hr (vertical and horiz relative distance) and angle of flow 
vr=bdepth-Dship+ynode; 

hr=((xnode-xchg) A 2+(znode-zchg) A 2) A 0.5; 

theta=atan((xnode-xchg)/(znode-zchg)); 

% 

% Free-field fluid velocities/accelerations (vertical and horizontal) 

% Horizontal and vertical free-field fluid velocities and accelerations 

dis=(vr A 2+hr A 2) A 0.5; 

uh=(el*hr/dis A 3)+(3*e2*vr*hr/dis A 5); 

u\ -(e 1 * vr/dis A 3 )+(3 *e2 * vr A 2/dis A 5 )-(e2/dis A 3 ); 

uhd=(e 1 d* hr/dis A 3)+(3 *e2d*hr*\T/dis A 5 )+. . . 

...vm*(3*el*hr*vr/dis A 5-(3*e2*hr/dis A 5)*(l-(5*vr A 2/dis A 2))); 

uvd=(e 1 d*vr/dis A 3+3 *e2d*\T A 2/dis A 5-e2d/dis A 3)-. . . 

...vm*(el/dis A 3-3*el*vr A 2/dis A 5+9*e2*vr/dis A 5-15*e2*vr A 3/dis A 7); 

% 

% Y and Z directed free-field fluid velocities and accelerations 
uv(n)=uv; 

uz(n)=uh*cos(theta); 

uvd(n)=u\d; 

uzd(n)=uhd*cos(theta); 

end 
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APPENDIX B 



"DRY" WHIPPING RESPONSE AT SPECIFIED NODES 
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acceleration (ft/sec2) 



Horizontal displacement vs. time for node 7 




Horizontal velocity vs. time for node 7 




Horizontal acceleration vs. time for node 7 




Figure B-l . Horizontal response at node 7, given initial displacments equal to a scaled 
horizontal mode shape vector for mode 1 flexure. 
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acceleration (ft/sec2) 



Horizontal displacement vs time for node 10 




Horizontal velocity vs. time for node 10 




Horizontal acceleration vs. time for node 10 




Figure B-2. Horizontal response at node 10, given initial displacments equal to a scaled 
horizontal mode shape vector for mode 1 flexure. 
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Horizontal displacement vs. time for node 20 




Horizontal velocity vs. time for node 20 




Horizontal acceleration vs. time for node 20 




Figure B-3. Horizontal response at node 20, given initial displacments equal to a scaled 
horizontal mode shape vector for mode 1 flexure. 
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Horizontal displacement vs. time for node 7 




Horizontal velocity vs. time for node 7 




Horizontal acceleration vs. time for node 7 




Figure B-4. Horizontal response at node 7, given initial displacments equal to a scaled 
horizontal mode shape vector for mode 2 flexure. 
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Horizontal displacement vs. time for node 10 




Horizontal velocity vs. time for node 10 




Horizontal acceleration vs. time for node 10 




Figure B-5 Horizontal response at node 10, given initial displacments equal to a scaled 
horizontal mode shape vector for mode 2 flexure. 



143 



Horizontal displacement vs. time for node 20 




Horizontal velocity vs. time for node 20 




Horizontal acceleration vs. time for node 20 




Figure B-6. Horizontal response at node 20, given initial displacments equal to a scaled 
horizontal mode shape vector for mode 2 flexure. 
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APPENDIX C 



"WET" WHIPPING RESPONSE AT SPECIFIED NODES 



145 



acceleration (ft/sec2) velocity (ft/sec) displacement (ft) 



Horizontal displacement vs. time for node 7 




Horizontal velocity vs. time for node 7 




Horizontal acceleration vs. time for node 7 




Figure C-l. Horizontal response at node 7, given initial displacments equal to a scaled 
horizontal mode shape vector for mode 1 flexure. 
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Horizontal displacement vs. time for node 10 




Horizontal velocity vs. time for node 10 




Horizontal acceleration vs. time for node 10 




Figure C-2. Horizontal response at node 10, given initial displacments equal to a scaled 
horizontal mode shape vector for mode 1 flexure. 
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acceleration (ft/sec2) 



Horizontal displacement vs. time for node 20 




Horizontal velocity vs. time for node 20 




Horizontal acceleration vs. time for node 20 




Figure C-3. Horizontal response at node 20, given initial displacments equal to a scaled 
horizontal mode shape vector for mode 1 flexure. 
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Horizontal displacement vs. time for node 7 




Horizontal velocity vs. time for node 7 




Horizontal acceleration vs. time for node 7 




Figure C-4 Horizontal response at node 7, given initial displacments equal to a scaled 
horizontal mode shape vector for mode 2 flexure. 
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Horizontal displacement vs. time for node 10 




Horizontal velocity vs. time for node 10 





0 0.5 1 1.5 2 2.5 3 

time (sec) 



Figure C-5. Horizontal response at node 10, given initial displacments equal to a scaled 
horizontal mode shape vector for mode 2 flexure. 
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Horizontal displacement vs. time for node 20 




Horizontal velocity vs. time for node 20 




Horizontal acceleration vs. time for node 20 




Figure C-6 Horizontal response at node 20, given initial displacments equal to a scaled 
horizontal mode shape vector for mode 2 flexure. 
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APPENDIX D 



WHIPPING RESPONSE FOR HULL SUBJECTED TO CHARGE DESIGNED TO 
EXCITE HORIZONTAL MODES 1 AND 2 FLEXURAL WHIPPING 
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Acceleration (ft/sec2) 



Horizontal displacement vs. time for node 7 




Horizontal velocity vs. time for node 7 




Horizontal acceleration vs. time for node 7 




Figure D-l . Horizontal response at node 7, for hull subject to bubble pulse designed to 
excite horizontal mode 1 flexural whipping. 
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Horizontal displacement vs. time for node 10 




Horizontal velocity vs. time for node 10 




Horizontal acceleration vs. time for node 10 




Figure D-2. Horizontal response at node 10, for hull subject to bubble pulse designed to 

excite horizontal mode 1 flexural whipping. 
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Horizontal displacement vs. time for node 20 




Horizontal velocity vs. time for node 20 




Horizontal acceleration vs. time for node 20 




Figure D-3. Horizontal response at node 20, for hull subject to bubble pulse designed to 

excite horizontal mode 1 flexural whipping. 
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Acceleration (ft/sec2) 



Horizontal displacement vs. time for node 7 




Horizontal velocity vs. time for node 7 




Horizontal acceleration vs. time for node 7 




Figure D-4 Horizontal response at node 7, for hull subject to bubble pulse designed to 
excite horizontal mode 2 flexural whipping. 
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Acceleration (ft/sec2) 



Horizontal displacement vs. time for node 10 




Horizontal velocity vs. time for node 10 




Horizontal acceleration vs. time for node 10 




Figure D-5. Horizontal response at node 10, for hull subject to bubble pulse designed to 

excite horizontal mode 2 flexural whipping. 
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Horizontal displacement vs. time for node 20 




Horizontal velocity vs. time for node 20 




Horizontal acceleration vs. time for node 20 




Figure D-6. Horizontal response at node 20, for hull subject to bubble pulse designed to 

excite horizontal mode 2 flexural whipping. 
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APPENDIX E 



WHIPPING RESPONSE FOR HULL SUBJECTED TO CHARGE DESIGNED TO 
EXCITE HORIZONTAL MODE 2 FLEXURAL WHIPPING AND INCLUDING 

INTERNAL MASS/ABSORBER 
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Horizontal displacement vs. time for node 7 




Horizontal velocity vs. time for node 7 




Horizontal acceleration vs. time for node 7 




Figure E- 1 . Horizontal response at node 7, for hull subject to bubble pulse designed to 
excite honzontal mode flexural whipping and including 200 LT "tuned" mass damper at 

node 7. 
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Acceleration (ft/sec2) 



Horizontal displacement vs. time for node 10 




Horizontal velocity vs. time for node 10 




Horizontal acceleration vs. time for node 10 




Figure E-2. Horizontal response at node 10, for hull subject to bubble pulse designed to 
excite horizontal mode flexural whipping and including 200 LT "tuned" mass damper at 

node 7. 
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Horizontal displacement vs. time for node 20 




Horizontal velocity vs. time for node 20 




Horizontal acceleration vs. time for node 20 




Figure E-3. Horizontal response at node 20, for hull subject to bubble pulse designed to 
excite horizontal mode flexural whipping and including 200 LT "tuned" mass damper at 

node 7. 
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Horizontal displacement vs. time for internal absorber 




Horizontal velocity vs. time for internal absorber 




Horizontal acceleration vs. time for internal absorber 




Figure E-4. Horizontal response of 200 LT "tuned" mass located at node 7 
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APPENDIX F 



NATURAL SLOSHING FREQUENCIES AND SLOSHING MODAL MASS FOR 
LATERAL SLOSHING OF RECTANGULAR TANKS 



164 



Weight of modal mass (Iton) Natural frequency (rad/sec) 



20 - 

15- 

10 - 

5- 

o- 

0 



2 3 4 5 6 

Lateral sloshing mode (n) 



60 

40 

20 

0 



0 



JS2- 



-S&L 



js&l 



2 3 4 

Lateral sloshing mode (n) 



5 



6 



Figure F-l . Natural sloshing frequency and sloshing modal mass for lateral sloshing of a 

(10'hx 10'wx lO'l) rectangular tank. 
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Figure F-2. Natural sloshing frequency and sloshing modal mass for lateral sloshing of a 

(20'h x 10'w x 20'1) rectangular tank. 
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APPENDIX G 



STRAIN AND VELOCITY HISTORIES (MEASURED AND CALCULATED) FOR 
THE "RED SNAPPER" MODEL TESTS 
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Stbd axial strain vs. time, element 6 




Stbd axial strain vs. time, element 9 




Time (sec) 



Figure G-l Starboard-side fiber strain histories for elements 6 and 9 for Test 2 of 
the "Red Snapper" model test (horizontal and vertical mode 1 whipping), and strains 
calculated using the computational algorithm. 
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Crown axial strain vs. time, element 6 




Keel axial strain vs. time, element 9 




Figure G-2. Crown and keel fiber strain histories for elements 6 and 9 (respectively) for 
Test 2 of the "Red Snapper" model test (horizontal and vertical mode 1 whipping), and 
strains calculated using the computational algorithm. 
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Velocity 



Horizontal velocity vs. time, node 10 




Horizontal velocity vs. time, node 20 




Time (sec) 



Figure G-3. Horizontal velocity histories for nodes 10 and 20 for Test 2 of 
the "Red Snapper" model test (horizontal and vertical mode 1 whipping), and velocities 
calculated using the computational algorithm. 
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Strain Strain 



Axial strain vs. time, element 6 




Axial strain vs. time, element 17 




Time (sec) 



Figure G-4. Starboard-side fiber strain histories for elements 6 and 17 for Test 3 of 
the "Red Snapper" model test (horizontal mode 2 whipping), and strains calculated using 

the computational algorithm. 
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Velocity Velocity 



Horizontal velocity vs. time, node 10 




Horizontal velocity vs. time, node 20 




Time (sec) 



Figure G-5. Horizontal velocity histories for nodes 10 and 20 for Test 3 of 
the "Red Snapper" model test (horizontal mode 2 whipping), and velocities calculated 

using the computational algorithm. 
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